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A rigorous derivation of the density functional in the Hohenberg-Kohn theory is 
presented. With no assumption regarding the magnitude of the electric coupling 
constant e"^ (or correlation), this work provides a firm basis for first-principles cal- 
culations. Using the auxiliary field method, in which need not be small, we show 
that the bosonic loop expansion of the exchange-correlation functional can be re- 
organized so as to be expressed entirely in terms of the Kohn-Sham single-particle 
orbitals and energies. The excitations of the many-particle system can be obtained 
within the same formalism. We also explicitly demonstrate at zero-temperature 
the single-particle limit, the weak-coupling limit of the energy functional, and its 
application to homogeneous electron gas. 

PACS numbers: 71.15.Mb 

I. INTRODUCTION 

At low energy scale, interactions among electrons largely determine the structure, 
phases, and stability of matter. Although this fact is well known, pragmatic first- 
priniciples/quantum- mechanical calculations to determine various properties of many- 
electron systems are often hindered by two factors. First, in most condensed matter systems, 
the typical interaction energy between electrons (the electric coupling constant divided 
by average electron-electron separation) is often larger than the typical kinetic/Fermi en- 
ergy of electrons. The results is that a perturbative expansion using as the expansion 
parameter may not be fruitful. This is particularly true for strongly correlated systems. 
Second, there is an exponential increase in the number of degrees of freedom as the number 
of electrons involved increases. When the number of electrons becomes large, according to 
Kohn,^ calculations based on constructing many-electron wave functions soon lose accuracy 
and will be stopped by an "exponential wall" . It is thus imperative to have a method that 
goes beyond the conventional perturbative scheme using as the expansion parameter and 
whose computational complexity does not grow exponentially with the number of electrons 
involved. 

In 1964, Hohenberg and Kohn^ proved a theorem stating that there exists a unique 
description of the ground state of a many-body system in terms of the expectation value of 
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the particle-density operator. This theorem started the development of the density functional 
theory (DFT), which offers a possibility of finding the ground state energy Eg by minimizing 
the energy functional Ey that depends on the charge density n only: 

Eg^mmEy[n\, (1) 

n 

with V representing the external one-particle potential of the system. The electronic density 
Ug, which minimizes the energy functional is the ground state electronic density. 

Hohenberg and Kohn showed that the energy functional Ey[n] can be decomposed into 

Ey [n] ^ J ^{r) n{r)+T [n] , (2) 

with T [n] being a universal functional independent of the external potential v. Mermin^ 
extended this theorem to finite temperature with Ey in (2) replaced by the grand potential, 
and J-[n] replaced by a different universal functional. The electron density ut, minimizing 
the grand potential functional, corresponds to the electron density at thermal equilibrium. 

To make practical use of the DFT, however, a recipe to compute the energy functional is 
needed. Kohn and Sham^ proposed a decomposition scheme, aiming to express the energy 
functional Ey [n] via an auxiliary, noninteracting system that yields a particle density identi- 
cal to that of the physical ground state. For a typical nonrelativistic many-fermion system, 
described by the Hamiltonian 



e 



the Kohn-Sham decomposition takes the form 



Ey [n] = To [n] + J Vion{x) n (x) dx - /iNe 



4j j'^^i.dy^E.M, (4) 

where the chemical potential ji is introduced to ensure / n(x) dtx. — N^, with being the 
number of electrons. Here Tq [n] is the kinetic energy of an auxiliary system of noninter- 
acting fermions that yields the electron density n (x) , and the density functional E^c [n] is 
the so-called exchange-correlation energy functional. Given E^c [n\ and provided that it is 
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differentiable, one may minimize the functional (4) to arrive at the famihar Kohn-Sham 
single-particle equations^ 

i=l 

All of the many-particle complexity is now completely hidden in the exchange-correlation 
energy functional. 

Although To[n] + -Ea:c[n] is universal,^ there exists no simple means thus far to obtain 
it. As a consequence, various ad hoc exchange-correlation density functionals have been 
suggested/needed to yield acceptable results in different settings®"^ when employing the 
Kohn-Sham scheme. Limitations of these approximate functionals have been discussed. ^'^'^^ 
Some of the failures while using ad hoc density functionals can be attributed to misuse of the 
density-functional theory. For example, it is sometimes neglected that the electron density 
n(x) achievable via introduction of a source potential must obey / n(x)(ix = Ne and does not 
cover the functional space {n(x) > 0}, a problem also known as the -u-representability.^^"^^ 
As pointed out in reference 15, neglecting these constraints may lead to conclusions^® that 
are not always valid. 

The main objective of the DFT is to describe a many-body system in terms of the 
expectation value of the particle density operator. In fact, the use of the expectation value 
of a suitable operator to describe a many-body system via Legendre transformation, first 
introduced into quantum field theory by Jona-Lasinio,^^ is known as the effective action 
formalism. As the temperature approaches zero, the effective potential becomes the ground 
state energy. This connection suggests that effective action formalism can be used to achieve 
the general goal of the DFT: describing a many-particle system in terms of the expectation 
value of the density operator. A number of publications^^"^^ showed that at zero temperature 
the effective action plus i^Ng is the ground state energy, hnking effective action to the DFT. 

Existing methods of expressing DFT via effective action formalism can be classified 
roughly into two categories: either (a) by introducing an auxiliary field or (b) by using 
a perturbative scheme assuming the electric coupling as a small parameter. The former 
category includes a method developed by Fukuda et al}^ and that developed by Polonyi 
and Sailer. The latter scheme is used by Valiev and Fernando^^ and by Fukuda et al}^. 
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The strengths of the auxihary field approach often come from the simphcity of the ef- 
fective action expression and from the fact that in principle each term already includes 
infinitely many Feynman diagrams. However, as pointed out by Fukuda et al.,^^ the aux- 
iliary field approach seems to lack a direct connection to the Kohn-Sham scheme. Valiev 
and Fernando^^ introduced an auxiliary field to compute the exchange-correlation energy. 
However, the source term they introduced is coupled to the auxiliary field instead of the 
electron density operator. Furthermore, as pointed out by the authors themselves, ^° an 
artificial decomposition of the auxiliary field into a sum of the Hartree potential, the ex- 
change correlation potential, and the remaining fiuctuations is needed to write down the 
exchange-correlation energy. 

There are two advantages when one uses electric coupling in the perturbative (diagram- 
matic) expansion without introducing an auxiliary field. First, under the effective action for- 
malism of this type, a direct connection to the Kohn-Sham scheme can be made.^^'^'^ Second, 
there exist other expansion-based developments that can be used to obtain TqIu] +Exc['n\, 
the universal density functional (UDF). For example, with increasingly complex incorpo- 
ration of KS orbitals and energies at each order of e^, Gorling and Levy^^ wrote Eg as a 
perturbation series. Employing the Luttinger-Ward^^ method that uses as the perturba- 
tive expansion parameter to calculate electron self-energy. Sham and Schliiter expressed E^c 
as a rather convoluted imphcit functional of electron density. ^^'^^ As shown by Tokatly and 
Pankratov,^^ these methods mentioned above can be expressed diagrammatically. However, 
the problem associated with based expansion remains. As described in reference 29, the 
expansion using is only good when is very small. Whether one can treat as a small 
parameter or not depends on the kinetic/Fermi energy of the electrons and the strength as 
well as the magnitude of variation of the single-particle potential involved. As a matter of 
fact, the success^° in employing the GW approximation^^ indicates that not treating as 
small may lead to results closer to experimental outcomes. 

In principle, the problem associated with assuming small can be tamed by summing 
an infinite subset of Feynman diagrams. However, as pointed out by Hedin,^^ it is nontrivial 
to devise a systematic resummation scheme where each new term is free from divergence 
even if ^ 1 and for which the sum of the new terms, each containing an infinitely 
many Feynman diagrams, accounts completely and non-redundantly for all conventional 
expansion diagrams. In this paper, without assuming small we develop an auxiliary field 
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method that makes a direct connection to the Kohn-Sham scheme and at the same time 
provides equivalently a systematic resummation scheme. 

A commonly used approximation for the density functional is the so-called local density 
approximation (LDA) in which the exchange-correlation energy is approximated by a linear 
functional Exc[n] ~ / drn{r) exc{n{r)), where exc{n) is a function of the local density (not a 
functional of the density profile) . See reference 1 for a nontechnical review. This approxima- 
tion ignores the nonlocal effect of the density profile, i.e., it assumes that 5Exc[n\/5n{Y) only 
depends on the value of n at r but not on the density n{r' ^ r) at locations other than r. To 
complement the LDA by incorporating nonlocal density dependence, Polonyi and Sailer^^ 
proposed the /-local approximation for the density functional, based on an idea very similar 
to the cluster expansion in statistical physics. Using this method to obtain explicit expres- 
sions for the approximate functional with Z > 3, however, becomes increasingly challenging 
due to the necessity of going through the coupling constant integration as required by the 
Hellmann-Feynman theorem.^^'^^ 

Another route to developing density functional is via the so-called optimal effective 
potential (OEP) methods. These methods typically start by introducing a priori an 
approximate, explicitly^^ orbital-dependent functional. (The approximate functional can be 
either Hartree-Fock or a more elaborated form.) The procedure then continues with a min- 
imization of the functional via varying single-particle KS orbitals and associated energies. 
Recently, the definition of OEP methods has been generahzed^^'^^ to include functionals 
dependent on either Green's function, the self-energy, or the KS potential. Since our ef- 
fective action based functional is based on self-consistently obtaining the KS potential, it 
falls exactly in the latter category. The generalized definition of OEP methods is probably 
becoming the standard definition now. A general characteristic of OEP methods is that the 
functional arguments -be they KS orbitals/energies. Green's functions, self energies, or KS 
potentials- are obtained via self-consistent procedure. Therefore, even with correction terms 
derived from pertrubative expansion, the self-consistency condition for OEP methods dis- 
tinguishes them from regular pertrubative methods. The important matter here is whether 
an OEP functional can be systematically improved and possibly be asymptotically exact or 
not. 

Based on effective action formalism, the OEP functional proposed here is asymptotically 
exact and can be shown to give rise to the desired UDF. Containing all the Mocal interaction 
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vertices, our method can provide equivalent approximate functionals for Z > 3 without 
going through the Hellmann-Feynman theorem. Since we focus on describing the proposed 
approach in a manner as self-contained as possible, we have included a non-negligible amount 
of standard materials available in existing literature/textbooks while keeping only a small 
portion of the existing literature that wc deem closely related to the present manuscript. 
Readers interested in gaining a broad background are referred to references 10 and 40-43 for 
reviews on the extensive body of literature in the DFT and related many-body approaches. 
Expert readers should note that new developments are mainly provided in sections IIIB- 
D. Although section IV and end of section V also contain some useful developments, they 
typically rcdcrivc/rc-cxprcss known results within our framework and/or provide contrast 
with existing methods. Section VI contains some insight of problems shared in post Hartree 
corrections. 

This paper is otherwise organized as follows. We first establish the notation in sec- 
tion II, followed by the development of the general formalism in section III. The purpose 
of subsection IIIG is to provide a computational recipe and to give some perspectives on 
computational complexity: no novelty is claimed here. In section IV, we discuss a num- 
ber of case studies: the emergence of the universal functional J^{n\ in Eq. (2) at arbitrary 
temperature, the behavior of the effective potential and the single electron limit at zero 
temperature, the screening effect, as well as the case of homogeneous electron gas. In sec- 
tion V, we then discuss the excitations of the system, and make comparisons with existing 
studies along this direction. An alternative formalism to obtain the effective action is then 
discussed in section VI. We conclude with the discussion and future directions section, in 
which we also provide some more relations/comparisons to other methods as well as some 
technical remarks. 

II. NOTATION 

Let us first define useful notation to lighten the exposition of the mathematical formulas. 
We define a three dimensional integral contraction by a dot 

a-h = J a(x) 6(x) 

where a and b may be single or composite fields. That is, with m > 1 and n > 1, a(x) and 
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6(x) may represent 

a(x) = ai(x) . . .a^(x) , 
and 6(x) = fei(x) . . . 6„(x) . 

Similarly, with a kernel M, we may define 

a-b-M = a-M-b = = dxciy M(x, y)a(x) b{y) . 

Note that in all expressions a is in front of b, and that is important when there are Grassmann 
variables involved. Evidently, one may generalize this notation to include higher order 
kernels. That is, one may have 



a-b-c-M = a-b-M-c = a-M-b-c = M-a-b-c ^JJJ dxdydz M(x, y, z) a(x) b{y) c(z) 
We define the four dimensional integral contraction by an open circle 

aob = J dx a{x) b{x) , 



with 



X = (r, x) , 
< T < ^ , 

r-/3 



and J dx — J dr J dx 



Again, a and b may be single or composite fields. That is, with m> 1 and n > 1, a{x) and 
b{x) may represent 

a{x) = ai{x) . . .am{x) , 
and b{x) — bi{x) . . . bn{x) . 

Similarly, we may define 

aoboM = aoMob = Moaob = J J dxdy M{x, y)a{x) b{y) , 

and 

Mottoboc = / / / dxdydz M{x, y, z) a{x) b{y) c{z) . 
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III. RELEVANT FORMULATION 

Consider the following generic fermionic Hamiltonian with s denoting the spins 



4e 

s.s' 



it(x)^;,(y)t/(x - y)4'(y)4(x)c^xdy, (7) 



+ 2 



^ ^ J ^li^) + '^ion(x, S) - - fXs^ V',(x) 

From this point on, we absorb — into i;ion(x, s). 

To lighten the notation, we will first ignore the spin degree of freedom but will comment 
on its effect when clarifications are needed. Let /3 be the inverse temperature. The partition 
function Z = Tr[e~'^^] contains all the information one needs. To probe the system in 
terms of the particle density, one often introduces a classical source term J(x) coupled to 
^^(x)^(x). The partition function now becomes a functional of the source J, and we write 



Z[J] e-^'^M = Tr 

It is easy to show that 



5W[J] 



Tr 



'^t(x)V;(x)e-^[^+'^-('^''^)] 



= (n(x))^ ^ nj(x) (9) 



5J(x) Z[J] 
Eq. (9) expresses n in terms of J, or more generally expresses (charge density of spin s) in 
terms of sources Jg' of all spins. Given fion(x) and C/(x — y) (for Coulomb interaction C/(x — 
y) = l^r^)' time-independent configuration of {J(x)} generates a time-independent 
charge density distribution {n(x)}. However, it is not guaranteed that every configuration 
of {r2(x)} is reachable by varying J. The functional variation on {n(x)} is thus limited to 
the subset of {n(x)} reachable by considering various {J(x)}. In this stationary case, the 
effective action r[n] is defined as the Legendre transformation of 

r[nj] = W[J] - J-nj , 
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where the subscript J indicates that the domain of r[n\ is the set of density profiles reachable 
by varying J, or the so-called v-representable^^"^^ densities. 

We now show the equivalence between the effective action and the energy functional Ey[n\ 
in (1). Since 



at zero temperature limit W[J] is simply the ground state energy corresponding to the 



while the electron density nj(x) is obtained by integrating all but one spatial variable of 
the ground state wave function corresponding to Hj. Evidently, when J — 0, r[n\\n=ng — 
l^[J]|j=o = Eg where Eg stands for the ground state energy corresponding to H and Ug 
represents the electron density at the physical (J = 0) ground state. When J ^ the 
corresponding electronic density n[J ^ 0] is different from Ug, and T[nj] represents the 
expectation value of the original Hamiltonian H, calculated using the ground state wave 
function corresponding to a different Hamiltonian, namely, H + J- {11)^11)). Since nj ^ Ug, 
r[n]|n=n[j] > r[n]|„=„g by the definition of the ground state. This means that r[n\ reaches 
its minimum at Ug, and various other electron density profiles n[J] producible by introducing 
different J form the domain of argument for T[n]. Thus T[n] has all the properties attributed 
to the energy functional Ey[n] in (1). Since the theorem of Hohenberg and Kohn states that 
this functional is unique, it must be equal to Ey[n]. As we will show in section IV B, r[n] 
can be decomposed exactly in the same manner as in (4). 

Within allowable configurations of {n(x)}, if one is able to invert the relation (9) to 
obtain, say, J[n\, then the explicit construction of the effective action becomes possible. In 
principle, this can be done via an inversion method. -"^^ Using this scheme, Valiev and Fer- 
nando^° proposed a perturbative expansion in terms of to express the exchange-correlation 
functional as a sum of an infinite number of Feynman diagrams and the diagrams' derivatives 
with respect to the Kohn-Sham potential. 

We approach this problem from two different routes, both involving the introduction 
of an auxiliary field. ^^'^^ As will be described in secton VIA, the second route does not 



Hamiltonian Hj = H + J-^tp'^tp) 



Hj ^ / dxV'^(x) -— V' + (z;ion(x) + J(x)) - V^(x) 




(10) 



10 



have an exact correspondence to the Kohn-Sham decomposition, but has the advantage 
that the correction terms may be obtained without further functional derivatives. The first 
route, as will be described later in this section, gives a recipe equivalent to the Kohn-Sham 
decomposition, together with a way to calculate the exchange-correlation functional in a 
self-consistent manner. 

This auxiliary field approach was pursued in an earlier publication,^* but there it was 
concluded that it seems infeasible to make a direct connection to the Kohn-Sham scheme. 
Using the inversion method, we show explicitly how the connection to the Kohn-Sham 
scheme can be made. One advantage of the auxiliary field method is that each Feynman 
diagram here corresponds to the sum of infinitely many Feynman diagrams in standard pcr- 
turbative field theory calculations^^ such as used in reference 20. Subsections III A through 
III F detail the proposed approach. Subsection III G lays out the computational procedure 
to give some perspectives on computational complexity. 

A. Path Integral 

To accommodate a time-dependent probe and to deal with excitations, we express Z as 
a path integral over Grassmann fields and we have 

^-m{J\ ^ z[J] = J D^^D^ exp{-5 [ip\^] - Jo(V^V)} , (11) 



with 



S [V^t, ^] = ^lJKG^K^|; + ^ (V^V) o Uo {^^^) , (12) 



2 

where ■0*'^^ denote Grassmann fields with ■0(^)(^,x) = — ■0*^^)(O, x), and 

dr - ^ + i;ion(x) - (x) {x\x') = ( dr - ^ + i^ion(x) - fi ] 6{x - x') 



2m J V 2m 

and 

U(x, x') = U(x -x')^5(t- t') C/(x - x'), 

5{x — x') = 5(r — r')5(x — x') . 

For a time-independent source, J{r,yi) — J(x) (i.e., Jo(t/;^t/;) — J dxJ{^)'4>\x)'^{x)). For 
time-dependent probes, J{x) becomes r-dependent, and Jo[%l)^%l)) = J dxJ{x)%l)^{x)%l){x). It 
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is straightforward to verify that 

5m [J]) 



{ip\x)ip{x))j = {fi{x))j = nj{x) . (13) 



5J{x) 

This quantity is important for later development. 

The quartic fermionic interaction in (12) can be disentangled via introducing an auxiliary 
real field with = Note that 

1 = Vd^ j D(t> e-i<^°^°'^ = VditC/ j D(t> e--2(^+y>^<'t>+y) , (14) 

for an arbitrary field Y, provided that Y{x) and Y{x') always commute. Since U{x,x') is 
diagonal in r, it suffices that they commute at equal Euclidean times. Let us set Y{x) — 
iip^{x)ip{x), which satisfies the equal time commutation requirement, and then multiply (11) 
by (14) to obtain 

Z[J]= J D(t)DiP^DiP exp{-5 [(t),i)\i)]] , (15) 

where 

S [0, = -\^HU) + + V'^°G'(;-,„-ioj)°V' , (16) 

with 

'^(0-i[/-loJ)(^'^') = (^9^-^+1^ion(x)-/X + i([/o0)^ + J(X)^(5(X-X') . (17) 

If we make a change of variable 0' = — iU~^oJ and then rename 0' by (j), we may rewrite 
(15-17) as 

Z[J] = e5-^°^"'°-^ j D(l)D7jj^Di; exp{-^j , (18) 

Sj [(f), V] = -^Tr ln{U) + Uo(f) + i(f)oJ + ipKG^hiP , (19) 
{x\G^'\x') = G^\x,x')= (^dr-^ + vU^)-i^ + i{Uo(l)),^6{x-x'). (20) 

Integrating over the Grassmann fields in (18), we obtain an effective theory in terms of 

= Z [ J] = e^^° j Dcf) exp {-/ [0] - iJo^} (21) 



where 



m = -^TVln(C/) + l0o[/o0 - TVln(G-i) . (22) 
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Let us introduce a new notation VF<^[J] via 



J D(j) exp {-/ [0] - iJo(j)} . (23) 



We describe later how to evaluate (23) using well-developed functional integral techniques. 
Evidently, we have 

[J] = [J] - ho U-KJ . (24) 

The expectation value, nj{x), of the density operator in the presence of a source term J 
is given by 

where 

{U-KJ), = J dyU-\x,y)J{y), 

and the expectation value of the auxiliary field is defined by 

..^ Sj^W^) /D0(#(x))e-^[^^-^°^ 

z^(x) = {^<|>(x))J = -jj^ = , (26) 

providing a relation between J and iif. 

Eq. (25) tells us that at the physical limit J — >■ 0, is the same as n. Since n is a real 
number, this implies that the expectation value of 0(x) is an imaginary number, which also 
implies that when viewed in the complex plane of 0(x), the saddle point of the integrand is 
located where the 0(x)s are imaginary numbers. 

Now let us write down the effective action. At finite temperature, the effective action is 
defined as the Legendre transform of /3H^[J]: 

r[n] = /3W [J] - ^i^M^l „ J = /3W4J] - -Jo U-KJ - noj . (27) 
Note that the functional derivative of r[n] with respect to n reads 



dT[n] 
5n 



'dif3W[J]) _ 



I- J- -J, (28) 



5J 

because 5{PW[J])/5J = n by Eq. (13). The effective action formalism requires one to express 
the probe J in terms of an expectation value of some sort, such as the electron density n, 
classical field i</7, or some equivalent quantity. Below, we will first calculate /31^^[J], and 
then use the system's electron density as the variable and make an explicit connection to 
the Kohn-Sham decomposition. 
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B. Evaluation of e-'^^'^I-^l via one-particle irreducible diagrams 

As shown by Jackiw,^^ it is possible to express VF<^[J] as a diagrammatic expansion con- 
taining only one-particle irreducible diagrams. The main idea is to shift the field by (/?, 
— >■ </? -I- 0, and note that J is a functional of </? via (26). We then rewrite 



where 



leading to 



Note that 



■]nZ,[J] = pW,[J] = - In 



/ 



D(f) e 



-I[(t>+<p]+I[<p]-iJo4) 



(29) 



(30) 



6J{x) 



leading to 



51 



S{i(p{y)) S{iip{y)) 

'<^'*-> - -J(y) . 



5J{x) 



Using an implicit method and replacing —J in (29) by the left-hand side (LHS) of (31), 
Jackiw^^ showed that is the sum of all connected one-particle-irreducible (IPI) vac- 

uum graphs governed by the action 

5I[p] 



+ + /[(^] +0O- 



6ip 



To evaluate the expression above, we first rewrite (20) as 

G~.U^,x) = G-\x,x')+iS{x - x')b{x) ^ G-\x,x') + V{x,x') , 



(32) 



with 



b = C/o(/) , 



and G^^{x,x') 
We may then write down 



dr- h t^ion(x) - /X + Uo{tip) 

2m 



6{x — x') 



(33) 
(34) 
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and 



Note also that 



In (G^) = In {G-') + Y: 



\k-l 



k=l 



k 



(35) 



[G^oY]^^^ = jdyG^{x,y)V{y,z) = J dyG^{x,y)5{y - z) {ib{y)) = G^{x,z) {ib{z)) 
Consequently, 

IVln(G'^y = IVln(G-i)+ j dx,G^{x,,x,){ih{x,)) 



-- I dxidx2G^(xi,X2)G^(x2,xi)(ib(xi))(ib(x2)) 



+ ^ J dxi. . .dxkG^{xk,Xi) . . .G^{xk-i,Xk){ib{xi)) . . .{ib{xk)) . (36) 



Therefore, 



with 



+ ^]+ I[ip] + <j) o^3^ = - ^6 oV-K b + y 7(^) [ip]o bi...obk 



fe=3 



p-1 ^U'^ -D , 
D{x, y) = G^{x, y) G^{y, x) , 



(37) 

(38) 
(39) 



and 



/^^^[(^]o bi...obk = ^ dxi... dxkG^{xk, xi) . . . G^{xk-i,Xk) [{ib{xi)) . . . {ib{xk))] . 

(40) 

As a side note, the quantity D{x, y) defined in (39) is also called the polarization associated 
with the Green's function G^, since it can be shown, by using a derivation identical to that 
leading to (71), that D{x,y) = —5G^{x,x)/5J{y) represents the reaction rate of density 
(given by n{x) = —G^{x,x)) due to the infiuence of the potential. According to Jackiw's 
results, Eq. (37) means that /3VF*[<^] is given by 



1 1 

PW^iip] = Trln{U) + -Tr In (p-^) - 

n=l 



,fe=3 



IPI, conn. 



(41) 



where the Trln(C/) term comes from the Jacobian of changing the variable from to 6 in 
(29), the angular bracket indicates the following average 
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and the subscript "IPI, conn." means to include only connected, one-particle-irreducible 
diagrams. In our context, a one-particle-irreducible diagrams refers to a diagram that cannot 
be separated into two by cutting a propagator line representing V. 
Substituting (22) and (41) into (30), we obtain 



1 °° 1 

+1tvi„(p-..c/)-E-L( 

n=l 



J^7«[(^]o6i...o6, 



k=3 



'IPI, conn. 



(43) 



Note that Fukuda et al}^ obtained an expression similar to (43) and used it to derive an 
effective action as a functional of i(f, which coincides with nj only at J = 0. 

We wish to keep nj as the functional variable. Let us first note that /3Vr[J] = /^^^^[J] — 
\JoU~^oJ . By employing the identity (25) 

nj — i(p — U~^oJ , 

we can now write /3VF[J] as 

/3W[J] = -Tt\n{G-')-^njoUonj + ^Tt\n(T)-KU^ 

oo " oo " " 

)lPI,conn., (44) 

n=l 

and 



,fc=3 



r[n] = -IVln(G-^) 



J + ^njo U 



on J + ^Trln ^V'^oU^ 



oo ^ 



n=l 



k=3 



)lPI, 



(45) 



For later convenience we introduce a parameter A to denote the order. Specifically, we write 



PW[J] = -Trln {G-') - ^njoUoUj + ^IVln (v-^U^ 



oo oo 

n=l ■ L k=3 

oo 



) IPI, conn. 



(46) 
(47) 



The bookkeeping parameter A will be set to 1 in the end. The exponent associated with 
the parameter A plays the role of the number of loops introduced. For example, the term 
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consists of one-loop contributions. The last part of (46) contains diagrams 
of two loops or higher. The explicit appearance of the njoUoUj term in (44) and (46) 
suggests the possibility of a connection to the Kohn-Sham decomposition of the DFT. 

C. Inversion Method by Loop Order 

We now come to the point of departure from typical treatments of auxiliary field approach 
and are ready to make a direct connection to the Kohn-Sham scheme. Let us first define a 
new free fermion propagator 



dr- h l^ion(x) - /i + Jo{x) 

2m 



6{x — x') , 



(48) 



where Jo is chosen such that the free fermion system has the same particle density as the 
physical system (where Coulomb interactions exist) with source potential J 



-Qo{x,x) = nj{x) . 



(49) 



The existence of Jo is scrutinized below. Prom the perspective of the Kohn-Sham decompo- 
sition, Eq. (49) corresponds to Eq. (6). In the presence of a source term J, Eq. (10) shows 
that it is equivalent to making Vion Vion + J- Comparing with Eq. (5), we see immediately 
if one were to choose 



Jo = J + U-nj + 



Sn 



(50) 



n=nj 



the requirement (49) can be fulfilled. We will therefore call the corresponding free fermion 
system the Kohn-Sham system. The presence of Jq, of course, depends crucially on the 
differentiability of -E2:cM, whose existence (not differentiability) was proven.^ 

To bring out the Kohn-Sham quantities (orbitals and energies) in our loop expansion, let 
us first use a variant of (25) 

uo(i(fi) — J + Uouj 

and replace Uo{i(p) by J -|- UoUj in the expression of propagator G^. Specifically, we write 
(34) as 



G-\x,x')^ 



dr - h Vion(x) - A* + J{x) + UoUj 

2m 



5(x - x') . 



(51) 



The critical step is to decompose the source J in a particular way. 



J = ( Jo - Uonj) + / = Jo + / . 



(52) 
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Therefore, from Eqs. (25), (34) and (52) we have 

G-\x, x') = Qo^x, x') + J'{x) 5{x - x') , (53) 

and the Kohn-Sham propagator Qq appears as we expected. As will be described in sec- 
tion HIE, Qq{x^x') can be expressed in terms of the Kohn-Sham quantities. Therefore, the 
idea is to expand around Qq provided that J' can also be expressed via Kohn-Sham 
quantities. Comparing (52) with (50), we find that formally speaking J'[n\ — — ^^l'^^^ ■ 

This source decomposition is introduced here for the first time in the auxiliary field 
approach. (A similar decomposition has been used^^'^° in the perturbative expansion in 
powers of e^.) 

The idea of inversion is to obtain J[n], that is, to find the corresponding J for each 
configuration of electron density n in the domain of V[n]. One then substitutes J[n] into 
r[n] — j3W[J[rii)[\ — J[n\on to express r[n] using the density profile n as the natural variable. 
For each given density profile n within the domain of F, Eq. (49) determines the correspond- 
ing Jq. The collection of such relations forms jo[n]. Similarly, if for every given n one can 
find the corresponding J', one obtains J'[n] and the goal of inversion is achieved. When 
evaluating r[n], one employs one density configuration at a time. That said, when we ex- 
pand H^[J[n]] = W^[Jo[n] + J'[n]] in powers of J' within the expression F[n] = — Jon, 
we will keep nj fixed, instead of treating it as a functional of Jq. 

We now examine the loop expansion of W[J\ carefully. Eq. (46) tells us that 

oo 

^W[J] = ^{Wo - W^o) + /9 5^ A' Wi[J + Uonj] , (54) 

1=0 

where I3{Wq - Wq) = -\njoUonj and I3Wq[J + UoUj] = -Trln(G'-^). Note that for any 
Wi term, its J dependence is through the propagator G^, which always has J -|- UoUj as 
the natural variable. Our reasoning earlier indicates that when we expand Wi[J -\-U on j] = 
Wi[Jq + UoTij + J'] = Wi[Jq + J'] in powers of J' within the expression F[n] = /3W^[J] — Jon, 
we may write the expansion in the following way (and forget about needing to keep nj fixed) 

Wi[Jo + J \ ^ Wi[Jo\-\ — — oj + o oj + ... . (55j 

OJq l\ OJqOJq 

With (55), we may express /3VF[J] as a double series 

PWIJ] = /3{Woo - Woo) ^ikJ''^' , (56) 

i,k 



18 



where each Wik involves the /c'th derivative of Wj. In particular, Wqo is given by (with 
nj ^ n hereafter) 

pWoQ = ^M^oo - \no Uon = -Tr HGo') - Uon , (57) 
and Wqi is given by 

PWoiiM = -Tr (^^^^^ ^-Jdzdygoiz,y)6{y - x)6{y - z) = -Qoix^x) = n , (58) 
in view of (49). 

Instead of looking at the expansion of Wi in powers of J' within the effective action 
expression, we now take a moment to look at the Wqq term in the double series expansion of 
^1^[J]. Consider the functional derivative of ^Wqo with respect to the source term Jq. Here, 
one is asking the response of ^Wqo with respect to change in Jq. Evidently, when Jq changes, 
its corresponding density n has to vary as well. Using the chain rule of differentiation and 
(58), we obtain 

— = — — — rio Uo^^ — no[l + (7o^^) — rio Uo^^ — n . (59) 

5Jo oJq SJq SJo SJo SJo 

This suggests that we define 

fo[n] = ^VFoo[^o] - Jo°n = -Trln(^?o"^) - ^no Uon - Jqou , (60) 

the Legendre transformation of /5 Woo [Jo], leading to 

5fo[n] 



Comparing (61) with (28), we find 

8{V\n\ - fo[n]) 



-Jo . (61) 



-J' . (62) 



5n 

The idea now is to develop a series for V\p\ led by Fofn]. Subtracting (60) from T\p\ — 
PW\J\ — Jon, we have 

V[n\ - f o[n] = ^W[J\ - ^#oo[Jo] - J' on , (63) 

in which the last two terms on the right-hand side (RHS) exactly cancel the terms in Woo 
and Woi contributing to ^W[J]. So the series for F — To is just (56) with those two terms 
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removed. Next we convert the double sum in (56) into a single sum by expanding J' as a 
series in A. We write 

oo 

J'N = 5]J,[r^]A^ (64) 
1=1 

where the precise expressions for Ji, J2, . . . are as yet undetermined since (64) is not a loop 

expansion. We substitute (64) formally into (63) and (56) to obtain a series 

00 

r[n]-fo[n] = ^r,[n] A^ (65) 
1=1 

in which each is defined explicitly in terms of the Jk, (3Wk<i[JQ], and their derivatives. 
Because Wqi is missing from (63), any occurrence of Jk is accompanied by at least one other 
factor Jk/ or else by an occurrence of some Wj>o, and hence by a power of A higher than the 
/c'th. In other words, the expression for Ti^i involves only Jk with k < I. We finally remove 
the indeterminacy in (64) by imposing (62) order by order in A, leading to (for / > 1) 

(-) 

Since r;>i involves only Jk<i, all the J; and F; can be found explicitly by applying (65) 
and (66) alternately. Evidently, it is the source decomposition (52) that allows us to obtain 
exact correspondence to the Kohn-Sham scheme. Below we will provide an explicit formula 
for Ti[n] in terms of VF;[Jo] and their functional derivatives. 

To obtain an explicit expression for Ti[n], let us substitute the LHS of (63) by the RHS 
of (65) and apply (56) as well as (64) to the RHS of (63). Then, by equating the coefficients 
associated with A' on both sides of (63), we obtain (for I > 1) 

SJn 



Fi [n\ = piWi [Jo\ + jj oJk 



k=l 

I ^ k,+Ml Sm^pWi.^,^^...^k^) [Jo]) 

^ — 5J0...SJ0 — ^67) 

m=2 fei,...,fe^>l " " 



For / = 1, we see that 

Fi[n] = /3Wi[Jo] = ^Trln (t)-^j^oU^ = Trln (v^'oU'^ . (68) 

We observe that 
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We then have 

J ^ JVM ^ 5.h 5{mi[M) 
^ 5n dn SJq 

Note that —6Jo/dn can be written as 

SJojx) ^ _ f 5n{y) \ ^ _ f P{lWMy' ^ ^^ro[n] 
5n{y) \5Jo{x)J \5Jo{x)5Jo{y) J 5n{x)5n{y) 



(69) 



where 

ToH = l3Wo[M - Joon = fo[n] - ^no Uon , (70) 
and the inverse is in the functional matrix sense. Using (49), we can evaluate Sn{x)/SJo{y) 

by 

= Qoix, y)Qo{y, x) = Dj_^j^{x, y) = Do{x, y) . (71) 

Therefore, 

Since —Qq[x, x) = n{x) represents the electron density of the KS system, we will call Dq{x, y) 
the polarization associated with the KS system. Note that if one were to approximate the 
effective action T[n] by T[n] = To[n] + Ti[n], the displayed equation above is the OEP 
equation for GW-OEP^^, while eq. (68) is the corresponding exchange-correlation functional. 
Once Ji is known, one can find r2 [n] 

dJQ Z OJqOJq 

z oJq oJq 
and J2 now can be computed via —5V2[n]/5n 

"^^ - ° SJo + ° SJoSJo + 2^° °^JMSJ^ ' ' ■ ^ ^ 
The explicit expression of J2 leads to ^^'^ on. It should now be clear how the strategy 
goes. We are assuming that the functionals {VTjJo]} ^-^d their derivatives with respect to 
Jo are known. This information can indeed be obtained by standard, albeit tedious, many- 
body perturbation method. One then uses Eq. (67) to express Fi in terms of the known 
functionals and with k < I — 1. Once F; is obtained, one can then use Eq. (66) to obtain 
Ji, which then facilitates the calculation of F^+i via (67) and so on. 
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We now take a moment to organize the terms of effective action r[n\. Prom (65) and 
(70), we know that 

oo oo 

r[n] = f oN + J2 r«N = ToN + -no Uon + J2 r^N ■ (74) 

1=1 1=1 
The first term on the RHS of (74) indeed corresponds to the effective action of the free KS 
system, the second term is exactly the Hartree energy. It is thus natural for us to define the 
last part, sum of Ti[n], as 

oo 

r,,[n] = Y.^i[n] . (75) 

1=1 

At the physical condition (i.e., when the source is absent), we should have 6T/6n = — J = 
0. Knowing that 

— = - Jo = - Jo + Uon , 

on 

we conclude that at the physical condition, 

— — !^ = Jo = Jo - [/ on . (76) 
On 

Note that the potential Jo together with the original non-interacting part of the Hamiltonian 
leads to the exact particle density of the interacting system. This implies that Jo contains 
both the Hartree term and the exchange-correlation potential. Since Jq = Jq + f/on^, with 
U onT being the Hartree term, the term Jq plays the role of exchange-correlation potential, 
as evidenced by (76). The quantum mechanical effects are completely contained in ra;c[n], 
defined in (75). 

Note that a different density profile other than that corresponding to the physical ground 
state can be induced by introducing a nonzero J. Let us again call the corresponding density 
profile nj. When J 7^ 0, we learn from (10) that it is equivalent to replacing v hj v + J. 
Eq. (5) then tells us that Jo must contain J, the Hartree term and the exchange-correlational 
potential as shown in (50). In this case, one writes Jq = J+{Jo — J)+Uonj. With Uonj being 
the Hartree term, {Jq — J) must represent the exchange-correlation potential corresponding 
to the configuration nj, and 

^J^^Jq-J. (77) 

This then leads to 

5n ' 
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what we expected when J 7^ 0. 

In terms of real computation, since the n dependence is through Jq, we may rewrite 
Eq. (76) as 

-r—0 — = Jq = Jq — u on 

On Jo 

or 

,^n„.l-'JE,j^^n„.(j„-u.n)-'-^^l^, (78) 

oJq dJo 
and the condition 6r[n]/Sn = is turned into Sr[Jo[n]]/6Jo = Since the effective action is 
a strictly convex function of the electron density n, there exists no local minima. One can 
therefore solve 5r[Jo[n]]/5n = by steepest descent. Effectively, we may define the direction 
of steepest descent k,{x) by 

^ jrmi = D„ o ( J„ - U,n) - llMmi (79) 

dJo{x) dJo 

and then update Jo{x) by Jo{x) — )■ Jo{x) + qn^x), with ^ > being the step size, till 
convergence is reached, that is, when k{x) 0. Note that eq. (78) is the standard OEP- 
equation for consistency. The iterative procedure described after (78) is largely identical to 
the Kuemmel-Perdew^^ procedure used for solving the exchange-only OEP equation. 

D. Diagrammatic Expansion of the Density Functional 

We now examine how one computes the effective action via diagrams. Prom Eqs. (65), 
(68) and (70), we have 

00 ^00 
r[n] = foN + J2 = ro[n] + -no Uon + J2 

i=l i=l 

11 °° 
= -Tr In (Qq^) - Joon + -no Uon + -Tr In (v^K +Y1 ' 

where 

p-i = C/-1 - Do (81) 
Do{x,y) = Go{x,y)Qo{y,x) . (82) 

Evidently, we need diagrammatic symbols for [/, and Vq. To get to higher-order terms 
ri>2 of the effective action, we see from Eqs. (71-73) and the text afterwards that it is 
necessary to incorporate into the diagrams the inverse density correlator D^^ and to evaluate 
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functional derivatives with respect to n (or Jq). We define tlie symbols for each line type 
below 



x' X U{x,x') x' ' =. a; D^^ 



x' . * 'X Qq{x,x') x'«^vvvwwn.x 1)0(3;, a;') . 

The smaller dots associated with the f/, -D,^^, and propagators are introduced to guide 
the eyes regarding the starting and ending points of these propagators. The propagator 
comes from contracting two h fields, see (33), and thus the bigger dots associated with V>q 
denote the coordinates/locations of those h fields. We will also use the bigger dots to indicate 
the space-time coordinate of a point of interest. 

To evaluate functional derivatives of Vri[Jo['^]] with respect to n (or Jq), we note from 
Eqs. (25), (34), (40) and (46) that the Jq dependence comes from Q<^{x^ y) and the functional 
derivatives associated with the formalism using Eqs. (67-79) necessarily require evaluations 
of SQo{x,x')/6Jo{y). Using a derivation similar to that in (71), we find that 

SGo{x,x') 

rjf. = -Go{x, y) Go{y, x ) . (83) 
(>My) 

The propagators D^^ and Vq also contain Qq and thus may be differentiated with respect 
to Jq. Note that Vq = (I — UoDq)^^oU and DQ{x,y) = Qo{x,y)Qo{y,x). Employing the 
identity 

— = -M-K—oM-^ , 

oJo oJo 

we let M = (I — UoDq) for the case of Vq and M = Dq for the case of Dq^ to obtain 

^ = vJ-^.V, , (84) 
and ^ = -D,K^-^oD,' . (85) 
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Eq. (83-85) may be expressed diagrammatically as follows 

X X 



6Dq\x,x') 



where the =F signs come from 



6 



6Qo(x,x') 



6Jo{y) SJoiy) 



X 
X 



SVo{x,x') 5 



6Jo{y) SJoiy) 



X' 
X 



6 



SJoiy) SJoiy) 



X 



X 



+ {y ) + iy 



' ^»(-') ^ 6 - 6 - 6 

z' z' z' 



SJoiy) 



SMy) 



When combined with the inverse density correlator, the graphs above yield 
6Qo{x,x') 



6n{z) 
6Vo{x, x') 
Sn{z) 

5Do\x,x') 
6n{z) 



-J dy DQ\z,y)go{x,y)go{y,x') , 

- J dydxidx2Df^^{z, y) Vq^x, xi) [^0(3^2, xi)go{xi, y)Go{y, X2)+ 
+Qo{xi,X2)Qo{x2, y)go{y, xi)] Vo{x2, x') , 
dydxidx2DQ^{z, y) Dq^{x, xi) [^0(3^2, xi)Qo{xi, y)Qo{y, ^2) + 
+^o(a^i, 3^2)^0(3^2, y)Go{y, xi)] Dq^{x2, x') , 
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which are shown diagrammatically below 

X 

5 



X 



5Qq{x,x') 

5n{z) 6n{z) 



X 
X 



X 



SVq{x, x') 
5n{z) 



6n{z) 



X 
X 



6Dq\x,x') 
6n{z) 



6n{z) 



+ 




+ 




The diagrammatic differential rules of 6/6Jq and 6/6n are needed not only for the cal- 
culations of higher-order terms Ti>2 of the effective action but also for the calculations of 
excitations, which we will discuss in section V. 

Before formally introducing the diagrammatic expansion, let us first set the convention 
that we will use. In general, each Feynman graph will carry with it a symmetry factor, 
which is inversely proportional to the number of ways to label this graph without changing 
the topology of the graph. The convention in quantum field theory usually leaves out 
the symmetry factor, as it may be deduced from the graph. To avoid enumeration of 
the symmetry factors needed, however, we will explicitly provide the symmetry factors for 
Feynman diagrams to be investigated later. 

Let us now look at the effective action (80) term by term. The Hartree term no U on/ 2 is 
expressed diagrammatically below 





P _1 

J- Hartree 2 

As remarked by Jackiw,^^ each term in the effective action expansion represents an infinite 
number of Feynman diagrams in regular perturbative field theoretic calculations. We use 
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the iTrln(Po of/) term as an explicit example 



2 ' '2 ' ' ^ 2n 

n=l 

where each term in the summation corresponds to a vacuum diagram, see Fig 1. 










Figure 1: Feynman diagrams corresponding to ^TilnCDQ^oU). The first term corresponds to 

— ^Tr{DooU), the second term corresponds to — -^Tr {Dqo U)'^ , the third term corresponds to 

— ^Tr{DooU)^, and the fourth term corresponds to — -^Tv {Dqo U)^ . In general, the diagram 
corresponding to — Tr (-Do° U)"' contains n bubbles strung hy n U propagators with the symmetry 
factor ^. 



If we pull together the lowest order diagrams in e^, we find the combination 

1 1 

-no Uon — -Tt{Dqo U) 

that corresponds to the Feynman diagrams in Fig. 2. 





(a) 



iP) 



Figure 2: The Feynman diagrams corresponding to the Hartree term (a) and the lowest order 
exchange term (b), the first graph from Fig. 1. The numerical factor associated with each diagram 
is shown explicitly. 



To compute F2, we need to first calculate Ji. Since ri[n] = /3I^i[Jo] = iTrln(I — DqoU) 



Mz) 



5Ti 



5n(z) 



-Tr 



[l-DooU)-\ 



SDn 



6n{z) 



-Tr 



Uo{I-DooU)-K-^^'' 



6n{z) 



-Tr 



Vno 



u 

6n{z) 



dy dx dx Dq ^{z, y)'Do{x, x')go{x, x')go{x' , y)Go{y, x) 



(86) 
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which can also be expressed diagrammatically as 

= ~ KD ■ 

From Eq. (72), we know that = f3W2[Jo] - iJioD^Kj^. We note from Eq. (46) and 

(47) that /3W^2[</o] corresponds to the diagrams in 



oo ^ ^ oo 



n\ ■ A 

71=1 L fe=3 



) 



IPI, conn. ) 



the last part of Eq. (46). There are two combinations of n and k that can give rise to A^. 
The first one is to have n = 1 and = 4, while the second one is to have n = 2 and k = 3. 
The first possibility generates two distinct graphs, while the second possibility generates 
three distinct diagrams. For n = 1 and k = 4, we have the following diagrams with their 
symmetry factors specified 

(a) (b) 





(!)!M)! 

1! 4 



1! 4 



For n = 2 and = 3, we have the following diagrams with their symmetry factors specified 




(c) 





3{{f 




-1)^ 


2! 


3 


3 


3(z)6 






2! 


3 


3 


3{{f 






2! 


3 


3 



The first diagram(a) among the three will be discarded since one can cut a Pq line and then 
separate it into two diagrams. Let us display below the diagram corresponding to JioDqoJi, 

JioDqoJi = { — Ji)oDqo( — Ji) = 






Therefore, the diagrammatic expression for r2[n] is given by 
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Figure 3: Diagrammatic expression for r2[ra]. 

Note that each Vq propagator hne contains the sum of infinitely many terms 

Vo = {U-^ - Do) = (I - UoDoy^ oU = U + UoDoo U + UoDqo UoDqo U+ ... (87) 




Since each U carries a factor of e^, the above expansion may be viewed as the expansion of 
T>Q with U being the leading order. Therefore, in the diagrams corresponding to r2 (Fig. 3), 
if one were to expand in e^, the first three diagrams are of order or higher, while the 
last two diagrams contains terms of order or higher only. 

It is now instructive to compare with the perturbative methods which use as the 
expansion parameter. Among those methods, the approach of Valiev and Fernando^" is 
the closest to ours. Let us pull out the diagrams contributing to order in rj<2 and 
compare to the results of reference 20. From ri[?T,], we see that the second diagram in Fig. 1 
corresponding to — |Tr [(Do°f^)^] wiH contribute to this order. The first three diagrams 
representing T2[n\ will contribute to the same order if we replace the V propagator by 
U . Thus, one obtains the diagrams shown in Fig. 4, which are identical to the results in 
reference 20. 




Figure 4: The Feynman diagrams of order U"^ (or e^) of the effective action T[n\. The correct 
symmetry factors are also provided. 



We use this example to illustrate the difference between our diagrammatic rules and 
those of reference 20. In order to obtain diagrams, the diagrammatic rules of reference 20 
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require the generation of all diagrams of /3iy[J], each of which is shown in Fig. 5. 




Figure 5: The Feynman diagrams of order U'^ (or e^) of /3VF[J]. The correct symmetry factors are 
also provided. 

Then one deletes diagrams that can be cut into two parts by cutting a Coulomb line u. 
This means that the last two diagrams above will be deleted from consideration. One then 
needs to find in the remaining diagrams the two-particle reducible (in the sense of fermion 
propagator) ones followed by an iterative operation to construct Dq^ lines. In this case, the 
only two-particle-reducible diagram is the third one and the iterative procedure generates 
exactly the only diagram containing a D^^ line in with T) replaced by U. Therefore, the 
diagrammatic rules of reference 20 require first one-particle- irreducibility of the Coulomb line 
followed by searching for diagrams that are two-particle-reducible (in fermion propagators 
sense) . 

For our case, when considering a diagram's reducibility, we only consider the Vq lines. Let 
us denote I^''^[Lp]obi . . . o bk by B^''\ and call it blob k. The one-particle-irreducible criterion 
simply means that when one joins any two blobs, say B^'^'^^ and B^''^\ one must make sure 
that at least two or more Vq lines are connecting B^''^^ and B^'^^^ due to contraction of b 
fields. Therefore, it is very easy to find and exclude one-particle-reducible diagrams in our 
method. For this particular example, only the diagrams surviving in the end are present in 
our formalism. Since the leading order of Vq, in terms of expansion of e"^, is U itself, our 
one-particle-irreducibility in T>q lines covers the one-particle-irreducibility of the Coulomb 
lines in reference 20. Although one may show that the rule^^ of inserting D^^ for two- 
particle-reducible (in terms of fermion propagators) diagrams still applies , it is not essential 
to have. However, one may choose to use this rule as a tool to ensure correct generation of 
all distinct diagrams. Equipped with the diagrammatic expansion rules, one may construct 
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r;>2 following the inversion method described in section IIIC. 



E. Evaluation of Qo using single particle orbitals 

To calculate Qo{x,x'), we define v(x) = Vion{^) — fJ' + Jo{^)- Since the evaluation of the 
Green's function is general, we use the symbol G{x, x') in place of Go{xi x') in the following 
derivation. 

Consider first a generic free fermion Hamiltonian, 



where 



(88) 



and the corresponding Green's function (with Z = Tre 



G{x,y) = (T^(a;)^t(^)) = ^-iTr [e-^^T(^(a;)^t(^)) 

The positive infinitesimal quantity r] is introduced to ensure that in the limit Tx — Ty, the 
time ordered product corresponds to normal ordering (represented by the second term at 
equal time). Note that < T^.Ty < (3, -f3 < - Ty < (3 and ■0(x) = e'^'^"'4^{-x)e~'^'^" and 
■0^(y) = e'^y^'il)^{y)e~'^y^ . When t^ — Ty < 0, one must have t^ — Ty + ^ > 0. Observe that 



G{tx -Ty<0) 



-Tr 
-Tr 
-Tr 



-G'(r,-T, + /3>0), 



(89) 



where the first equality results from Ti^AB) — TiiBA), while the final equality results from 
the definition of the Green's function with positive time argument T^^^ — Ty > 0. Similarly, 
one may easily show that G{tx — Ty > 0) = —G{tx — Ty — [5 < 0). Therefore, the Green's 
function is antiperiodic in the imaginary time r. 

Letting M{x, y) = (drj{x - y) + h{x, y)5(r^ - Ty)^ = (dr, - S + "^W) S{x - y), we 
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find 



8m my) " 

S 5 Jp[^t^^]e-V't°Mov-+^-ov-+^toi; 



S 



JoM-'^oS, 



M-\x,y). 



This implies that / dx'M{x,x')G{x',y) — (^dr^ + /ix^ G{x,y) — S{x — y) with /ix being a 
one particle first quantized Hamiltonian /ix = + '^(^)- Note that 5{x — y) — 5(x — 

y)8{.T~x — T~y) and the latter delta function in time is defined via g{Tx)5(Tx — ry)dTx = gijy) 
for any antiperiodic function 5'(r). With this understanding, one may express 5{x — y) in 
the following way to obtain G{x,y). Observing that 

{dr, + h^)G{x,y) = {x\y) = ^{x\un,a){un,a\y) 

-iojjTj-Ty) 



OJn,Ce 



one sees that the action of {dj-^ + h^x) may be compensated, leading to 



G(x,t/) = 5]0«(x)C(y) 



-y 



where 



(90) 



{x\uJn,C() 



^a{^)- 



7r(2n + 1) 



and /ix0a(x) = 
Eq. (93) imphes that 



2m 



+ i;(x) 



6a (x) = {Sa - lj)(t)a{y^) 



-7, ^'^ion(x) + Jo(x) 

2m 



0a (x) = £a0a(x) . 



(91) 
(92) 

(93) 



(94) 



Only frequencies of the type is included in the expansion to ensure the antiperiodicity 

of the fermionic Green's function. To proceed further in (90), one may sum the frequency 
by introducing a function —^/{e^'^ + 1) which has poles at a; = »7r(2n+i) residue 1. 
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Evidently, poles for the function —^jieP'^ + 1) occur whenever e^^ + 1 = 0. To investigate 
the strength of each pole, let's rewrite e^^ + 1 in the limit when a; — >■ i7r(2n + as 

(. -.X n, i7r(2n + l)A , ( ^, i7r(2n + l), 
1 + exp i7r(2n + 1) + ^(u; ^— '-) = 1 - exp ^— ^) 



= a; 



i7r(2n + 1) 



/3 



ijj 



i7r{2n + T 



Therefore, —/3/{e^'^ + 1) indeed has residue strength 1 at each of the allowable frequencies. 

To evaluate 4 , . "["^ — — when < t„, one integrates over a circular contour (with 
radius |a;| — )■ oo) on the complex w-plane 



-duj . 



Because the line integral along the infinite circle vanishes, the sum of residues must vanish, 
meaning that 

0. 



Lp-(^a-M)(rx-rj,). ^ 



-iujn + ea- u e^^^"-'^') + 1 

Thus, when Tx < Ty (because when Tx = Ty, it must agree with r^; — — )■ 0~) one finds 

G{x,y) = J]0„(x)0;(y)e-(^«-'')(-^--^)(-n«) . 

a 

On the other hand, when > Ty, one considers 



The residue sum then turns into 



E 



Q—i0Jn{Tx—Ty) 

-iuJn + ea- 1^ 



g-/3(£a-/i) _|_ \ 



. 



Thus, when > Ty one obtains 



Therefore, we have 



-{eoc- IJ.){Tx-Ty) 



'-Ua) if Tx < Ty 
1 - ria) if Tx > Ty 



(95) 
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where Ua — l/(e^^^"~^^ + 1). Note that in the expression involving £q,, it is always minus 
the chemical potential /i. 

To evaluate Qo{x,y), we need to solve the eingensystem (93). It requires evaluations of 
the RHS of Eq. (79) and self-consistency is reached when k{x) — > 0. Of course, one cannot 
evaluate all the terms and must truncate the series on the RHS of Eq. (79) at some stage. 
This implies that the density profile obtained through the self-consistency condition in this 
manner depends on the number of terms one includes on the RHS of Eq. (79). Nevertheless, 
the self-consistent solution obtained when keeping k terms on the RHS of Eq. (79) can be 
used as the starting point when one wishes to include A; -|- 1 terms on the RHS of Eq. (79). 

F. Functional derivative of Qo{x,x') and / dry Qo{x,y)Qo{y,x') 

Especially when Jq is time-independent, we need to evaluate 



SQo{x,x') 

SMy) 



SMy) 



Qq{z' , x') dzdz' 

— j Goix, z) 5{z — z')6{y — z) Qo{z',x') dzdz' 

- / So{x,y)go{y:x')dry, (96) 
Jo 



where the central expression So{x, y) Qoiy, x') dry may be re-expressed by single-particle 
orbitals as will be shown below. 
From Eq. (90), we see that 

1 „—iu)„/{Ty—T^i—ri) 

and G{y,x') = -J^ ,^ , , My^i^') , 

n',p 

and therefore 

■t-^ ( 1 g-ia;„(Ta,-T^/-27j) \ 

= E E H.„..„-.)(-^„ + J my)W(^) 

Let us now concentrate on the portion inside the parentheses. Assuming that there is no 
energy degeneracy, we rewrite the denominator of this fetor when a ^ p 

1 1/1 1 \ 



{-iUn + £a - lj){-iuJn + Ep - /d) Ep - Ea \ -lOJn + Ea - /jL -lOJn + Ep - H 



34 



and when a — p 



d 



Therefore, we need to evaluate ^ — — ■ 

We distinguish the two cases: < Tx' and r^; > Tx'- When Tx < Tx', we consider the 
following integral over the infinite circle 

-CO + £a- l^e^"^ + 1 2m ~ ^ -icOn + £a - A* e^(="-^) + 1 
Since the integral over the infinite circle vanishes, we have 

g-(£a-M)(T-x-r^/) 



= — e 



-(£a-M)(-ra:-v) 



(97) 



To evaluate the expression 4 . ."l — we consider 



E 



dea ^ + ea- fJ. 
Therefore, when Tx < Tx', 



e-(e„-M)(r.-v) [-(3n^^i _ _ ^rx - Tx>)n, 



^0 

(98) 



+ ^</.„(x)0;(y)0,(y)0;(x') 



where n„(y) = 0;(y)0«(y). 

On the other hand, when Tx > Tx', we consider the integral 



gU;(Tx-r^/) _p g-jw„(Ta,-r^/) ^g-(£«-M)(Ta;-T^/) 

— -)■ 



a; + - /X e/^'^ + 1 27ri ^ + - p e-/^(^«-'*) + 1 ' 
and obtain (since the integral over the infinite circle vanishes) 



1 „ Q-i^^n{Tx-T^l) 



3-(£a-/i)(rx-T^/) 



,-(£a-/i)(T^-T^,)/l _ 



(1 - n„) . 



(99) 



Similarly, 
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Therefore, when > t^', 



dry G{x, y) G{y, x') = J] 0«(x)n«(y)0;(x')e-(^"-'^)(-^--^') [-/3n,(l - n<,) + (r, - t,0(1 - n,)] 

Up e n« ^^^^^ 



We may now write down the full expression for SQQ{x,x')/SJQ{y) as 



SQo{x,x') 

SJoiy) 



- I Qo{x,y)Qo{y,x')dTy 
Jo 

na(y) C(x') e-^^"-'')^^^-^^') - n«) + (r, - r^On.] 



-Y.M^'^^^MMyWpi^') 

ay^p 

In the absence of time-dependence, we have 

1 S{f3W[Jo]) 



g-(£p-M)(Tx-v)^^ _ g-(£a-M)(-rx-T^/)^^ 



(101) 



n X 



/3 5Jo(x) 
Therefore, using Eq. (101) we have 

5n(x) 



= Jdzdygo{z,y)S{y - x)6{y - z) = -goix,x) . 



^p 
^p 



+ X]^a(x)na(y) [/3na(l-n«)] 



Note that the expression ^11^(1 — fia) vanishes as (3 — )■ oo, because 11^(1 — Ua) decays 
exponentially with /3 when /j, ^ That is, at the zero temperature limit, one has 



^p 



L'o(x,y) E<^,(x)0;(y)0p(y)0;(x) 
as long as /X is not equal to any eigenenergy of the orbital. 



G. The Computational Procedure 

The recipe for computation goes as follows. Starting with a reasonably guessed Jo(x), 
one first obtains single particle wave functions 0q and energies through (94). Note that 
the occupation number of state a is given by 

1 



g/3(£a-/i) _j_ Y 
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and the chemical potential is chosen such that 

gs^na = Ne , 

a 

where Qs denotes the spin degeneracy {qs = 2 for spin 1/2 electrons). 

Through eqs. (95), (82), and (81), one constructs respectively the Qq, Dq^, and Vq 
propagators. Combined with their differentiation rules (83-85) with respect to Jq, these 
propagators are used to compute 5V /5Jq. 

As the third step, one obtains a new estimate of Jo given by Jq — <;5V/5Jq with > 
being the step size. See eq. (79) and the text nearby for details. 

Finally, one starts again with the new Jq and goes through the other steps, iteratively 
until the convergence condition 5V/5Jq = is reached. This simultaneously determines Jq 
(sum of the Hartree potential and the KS potential), the ground state charge density and 
the ground state energy, as well as the KS orbitals and energies. 

IV. CASE STUDIES 
A. The Universal Functional T[n] at Arbitrary Temperature 

There are vast discussions-'^^'^^'^-^ on the equivalence between Ey[n\ in (1-2) and 
lim^^oo ^r[n]. However, not much attention was given to the emergence of the univer- 
sal functional J^\n\ (at any given temperature) resulting from the effective action formalism. 
We address here for the first time how J^[n\ arises naturally from our effective action formu- 
lation. 

Fukuda et al}^ showed that it is possible to eliminate the v dependence in the functional 
to obtain (when translated into our terms) 

They interpret i(p as some sort of electron density since it coincides with the real electron 
density when J = 0. The appearance of the quadratic term in v, however, disagrees with 
(2) where it is clearly stated that in addition to a term that is linear in v, the remainder 
should be -u-independent. We settle this discrepancy below by showing that the appearance 
of the quadratic v dependence is due to the fact that reference 18 did not use the system's 
electron density as the natural variable. Once the system's electron density is used as the 
natural variable, the universal functional T[n\ emerges naturally. 
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Note that (10) indicates that one may view fion(x) as a nonvanishing source term. In 
the change of variable — >■ + iU~^oJ right after (17), if we make (f) ^ (f) -\- iU~^o(J + v) 
instead, we see that v is then bonded with J from that point on. That said, we may view 
W[J] = W'[J + v]. That is, the generating functional in the presence of the one-body 
potential v with source J is equivalent to the generating functional of a system without a 
one-body potential but with source J + v. Following the algebra in Eqs. (21-46), one sees 
that 



where 



/3W[J] = /3iy;[J + v]- -{J + v)oU-\J + v) , 



pW'^[J + v] = -(^oC/o(^-Trln(G'-^)+i(J + t;)o(^ 



(102) 



1 °° 1 

n=l 



5;;/W[(^]o6i...o6fe 



,fc=3 



'IPI, conn. 



5r - — - A* i{Uoip)x ) S{x - x') . 



Upon using the following variant of (25) 

n — icp — U~^o{J + v) , 

the quadratic term in J + v in (102) gets absorbed into the density of the electron and one 
arrives at the following variant of (44) 

PW[J] = PW'[J + v] = -Trln{G^^) --noUon+-Trln('D-Ku) 

2 2 V / 



oo ^ 



n=l 



5^/W[99]o6i...o6fe 



,fe=3 



IPI, conn. 



Note that here n = nj = n'j_^_^ represents the electron density of the system. 
The effective action r[n] = — Jon can thus be rewritten as 

r[n] = I3W'[J + v]- Jon = von + I3W'[J + v]-(J + v)on , 

where the first term on the RHS represents P J n(x)t;(x)(ix, and the last two terms represent 
the effective action, r'[n] — f5W'[J + v] — {J + v)on, of a system in the absence of one-body 
potential. We thus identify T[n] as 

J^[n] = ^ {PW'[J + v]-{J + v)on} . 

This also serves as an alternative exposition of what Mermin proved.^ 
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B. Effective Potential near Zero Temperature 

The effective potential divided by ^ is the ground state energy plus jiN^ in the T — > 
limit. Below, we will show how the Hartree-Fock terms appear in this formulation. 

Equations (47-73) provide a systematic expansion for calculating the effective potential. 
In the expression (60), the term — Trln(^(^^) is equivalent to — ln[det(^Q"^)]. There are 
many ways to obtain det(^Q"^): one may obtain the results directly through the definition 
of 6"'^'^°''^°] or one may express e~^^''t'^°l as a path integral to obtain the determinant in 
discrete time. We shall denote — Trln(^(^^) by /3Wo[Jo] with 

= j P[^t^^]e-Vto0„-iov, ^ det(^?o"') • 
Direct evaluation using the trace definition leads to 

det(^o"') = n (1 + e-^^'--'^^) , 

a 

and consequently 

a a 

Note that the energy is measured with respect to the chemical potential ji. Therefore, in 
the limit of ^ — > oo, we have 




(103) 



with £a < A* for 1 < q; < A^e- 

Using Eqs. (80) and (103), we obtain the low temperature limit of the effective action 




(104) 



The first term in (104) can be expressed in a different way if we multiply both sides of 
Eq. (93) by 0q(x), sum a over the lowest states, and integrate over x. Upon doing this. 
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we obtain 



h-Uionlx) - //+ Jo(x) 

2m 



J V e J ' e /• 

a=l a=l 

= To[n] - jiNe + y ["^ionlx) + Jo(x)] n(x) , 



(x) 



where 



and 



a=l 



To[n] ^ 5^ / dx0;(x) 



2m 



0a (x) . 



Therefore, the first two terms in (104) give us 



^(£a - 1^) - Jo-n = To[n] - iiN^ + / cix'Uion(x) n(x) . 



Evidently, the third term in (104) is nothing but the Hartree energy 



1 TT 1 

-n- U-n = - 

2 2 



(ix(iyn(x)-j |'T'(y) ■ 



Therefore, we have 

^lirn^ ^-^rfn]^ = To[n] + J i;ion(x) n(x) o?x - //A^e + ^ jj c^xdy n(x)|— ^— |n(y) 



+ hm ^Tr In (I - Dqo U) + hm 

/J->oo 2p /3->-oo 



^ oo 



1=2 



(105) 



with the last two terms combined to form the exchange-correlation energy functional when 
compared with the Kohn-Sham decomposition (4). 

If Dqo U (or e^) may be treated as a small quantity, one may expand 2^Tr ln(I — Dqo U) 

as 

(106) 



n=l 



with the leading term (in the static limit) 



-1 



— j dxdyDo{x,y)u{y,x) = — 



dxdy 



n(x, y)n(y,x) 



|x-y| 
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where n(x,y) = Z)m=i ^«(x)0* (y) = -^o(x, r; y, r). Consequently, when DqoU (or e^) can 
be treated as a small quantity, we may write the leading terms of the effective potential as 

-'^ '' n(x)n(y) -n(x,y)n(y,x) 



1 /■ f 1 

lim -^r[n] = Tofn] + / Vion(x)n(x)dx - //A^e + / c^xdy- 

^^•oo jj J 2 J 

CO ^ ^ OO 

[(B„„ari + -.^r.i, 



|x-y| 



+ lim 

/3->-oo 



^ n=2 ^ j=2 



(107) 



where the fourth term is nothing but the Hartree-Fock term. The higher-order terms inside 
the square brackets encode the remaining contributions of exchange and correlation. In the 
case when the Coulomb interaction is strong, one may choose not to use the expansion in 
Eq. (106), but to use (105) for the zero temperature limit and (80) for finite temperature. 

C. Single Electron at Zero Temperature 

When the system contains only one electron and when the energy difference between the 
first excited state and the ground state is much larger than ksT, there will be no particle- 
hole pairs. Consequently, there should be no exchange-correlation energy as well as the 
Hartree energy. When one employs empirical density functionals, this feature is unlikely 
to be preserved - an issue known as the self-interaction problem. It is customary to define 
the exchange-correlation energy E^c as the sum of the Fock exchange energy and the 
correlation energy E^. Since it is easy to show that the exchange energy E^ exactly cancels 
the Hartree energy in the case of one electron, one easily concludes that Ec = for the one 
electron case. This argument has been used, for example, by Perdew and Zunger."^^ However, 
from the diagrammatic expansion point of view, the Hartree term and the Fock exchange 
term correspond only to the first order (in terms of e^) diagrams. The cancellation of the 
first order terms does not imply that the higher order diagrams, making up Ec, will give no 
contribution. That is, although the fact that Ec — for one electron system can be argued, 
a formal diagrammatic exposition incorporating all higher orders is needed to justify the 
asymptotic exactness of the proposed functional. 

The purpose of this section is to illustrate how E^. = can be derived formally within our 
formalism. It should be noted, however, that the self-interaction will remain if one truncate 
the series in eq. (80). While the exchange-only functional will have no self- interaction prob- 
lem, as we will show below it is not because the exchange-only functional is an approach 
with more correct physics, but because the simplification it makes is equivalent to completely 
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ignoring the correlation energy. 

When Nf, — 1 and when T — > 0, we find from Eq. (95) that the Green's function Qo{x, y) 
takes the following form 

' (/)i(x)0t(y)e-(^i-'^)(-=^---)(-l) if < Ty 
= ^ E 0.(x)C(y)e-(--)^--) if > ' ^'''^ 

a>2 

where Si < ji < ea>2 and 0i(p)(x) is the ground (p^^) state wave function of the single 
particle Hamiltonian 

h{x) = -— + l'ion(x) + Jo(x) - H 

with eigenvalue £i(p) — /x. 

The disappearance of the Hartree energy and the exchange-correlation energy is best seen 
by grouping terms with the same number of Coulomb lines U. We will explicitly show the 
first few calculations followed by a sketch of the general proof. 

Let us first show that the first term (n = 1) on the RHS of (106) cancels the Hartree 
term exactly. When there is only one electron, the Hartree term becomes 

f d^dy'^^'^Hy^ - /• (x)0t(x)0i(y)0t(y) 



J dxdy— 



|x-y| z J |x-y| 
The Fock term (n = 1 term of order U in (106)), when there is only one electron, reads. 

The cancellation between the Fock term and the Hartree term is thus apparent. If one 
were to approximate the exchange-correlation functional by the Fock exchange functional, 

correlation energy can never be accounted for. That is, this approximation coincidently 
leads to the expected result = at the single electron limit simply because it never takes 
Ec into account. 

In terms of diagrammatic expression, the Hartree term is given by diagram (a) of Fig. 2, 
while the lowest order exchange term {n — 1 term of (106) ) corresponds to diagram (b) 
of Fig. 2. That is, the order U diagrams in ^T[n] are identical to the order U diagrams 

in Vr[Jo] at zero temperature. The perfect cancellation between the Hartree term and the 
lowest order exchange term implies that the sum of these two terms remains zero when one 
makes a derivative with respect to either Jq or n. Diagrammatically speaking, this means 
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that when the system contains only one electron, we have 
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SJo _ ^ 


SJo 


1 


(50 


± 


2 


Sn 


^ 5n 








(109) 



Terms of the next order in U are contained in rj<2. See Fig. 4 for a diagrammatic 
expression of Ti<2 of order U"^ (or of order e^). We show below that in Fig. 4 the first two 
graphs cancel each other and the last two diagrams also cancel each other. For illustration, 
we will work out the explicit cancellation of the first two graphs by labelling the space-time 
points. The cancellation of the the last two graphs require additional operations which we 
will turn to later. The first two graphs in Fig. 4 with the space-time points labelled appear 
to be 

yi yi z/2 





^ ^ TT^ ^ Go{xi,yi){Qo{x2,y2) 

drJ,Ty dxidyi— 

-L-L 4xi-X2 yi-y2 



i=l 



>^[Go{yi,x2)Qo{y2,xi) - Go{yi,xi)go{y2,x2)] . 



(110) 



We have used to denote the time associated with both xi and X2, while denoting by Ty 
the time associated with both yi and y2. When Ty < r^., the quantity inside the square 
brackets of Eq. (110) vanishes because 

GoiyuXj) oc (f)i{y,)(j)l{yij) , 

and thus ^o(l/i, 2:2)^0(1/2, a;i) = Qo{yi,Xi)Qo{y2,X2)- When Ty > r^., we simply swap the 
dummy variable yi and ?/2 in the second graph above to arrive at 
2 



(iTarfr^]^ d:>^idyi ^^°jf^'^'^}^^^^^^'^'^y^ [Qo{xi,yi)Qo{x2, yi) - Qo{xi,y2)Qo{x2,yi)] ■ 



i=l 



4|xi - X2||yi - y2| 



And again the quantity inside the square brackets vanishes due to the same reason as before. 
As for the cancellation of the last two graphs in Fig. 4, we first use the bottom portion of 
(109) to obtain 



5n 







+ 
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and then 





or 



2 V^yv^ 2 








(111) 



when there is only one electron in the system. 

Equation (111) means that at the single electron limit, the last two graphs in Fig. 4 are 
equivalent to the last three graphs in Fig. 5. In other words, at the single electron limit and 
at zero temperature, the set of order U'^ diagrams in T[n\ is equivalent to the set of order 
f/^ diagrams in /3iy[Jo]. We shall pause at this point and elucidate the general situation by 
looking at these cancellations via Hugenholtz diagrams. 

With a two-body interaction term such as the Coulomb interaction, typical Feynman 
diagrams treat the direct (Coulomb) and exchange matrix element separately. It is not 
surprising that one may simplify the calculation by combining the direct and exchange 
matrix elements into a single antisymmetrized matrix element. The basic idea is to combine 
the following two scenarios into one 

7 5 7 6 1 S 






a pa pa p 

(jdlvlap) — {'y6\v\pa) = {■j6\v\ap} , 



;ii2) 



where 



and 



{-fS\v\ap)= / dT^dTydyidy(j)*Jx)(j)*g{y)v{x,y)(f)a{:>i.)(j)p{y) 



{-f6\v\ap} = J dr^c/rydxrfy0*(x)05(y)t;(x,?/) [0„(x)0p(y) -0p(x)0c,(y)] , 

with 0a (x) being the single-particle wave function described earlier in Eq. (93). The re- 
sulting diagrams with bullet dots as the new vertices are called Hugenholtz diagrams. The 
appearance of those vertex matrix elements comes from the following. When one evaluates 
a Feynman diagram, a propagator Qo{x, x') going from x' to x connects two vertices located 
at x' and x respectively. From Eq. (95), we know that 



Go{x,x') 



[-Ua) if < T^' 
[1 - ria) if > T^' 
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Upon integration of the space-time coordinates, we see that 0q,(x) will be integrated with 
vertex x, whither the propagator leads, while 0* (x') will be integrated with vertex x', whence 
the propagator originates. Therefore, when evaluating a Feynman diagram, one may asso- 
ciate each propagator going from time r' to time r with a state index a with 

Q^{a, t-t') = e-(--'^)(---') [(1 - n„)^(r - r' - r/) - nJ{T' - t + r])] . 

Each vertex will contribute a numerical factor equivalent to its vertex matrix element. Note 
that each vertex carries a time index. At zero temperature, for a vertex with time index r, 
if the two incoming lines originate from vertices with times larger than or equal to r, the 
vertex matrix element associated with r becomes zero in the single electron limit. This is 
because both incoming lines each carry only the a = 1 index. Upon antisymmetrization, 
the Hugenholtz vertex matrix element becomes zero. 

For an arbitrary Hugenholtz diagram with n vertices, one may always name their time 
index such that ti < r2 < rs < ■ • • < r„. In this case, the vertex associated with ti gives 
zero matrix element since its two incoming lines must come from two other time indices 
that are larger than or equal to ti. Consequently, each Hugenholtz diagram yields value 
zero at zero temperature when there is only one electron present. As a matter of fact, the 
order U diagrams of W^[Jo] are represented by the Hugenholtz diagram iQO- C)n the other 
hand, the first two graphs of Fig. 4 (order f/^ terms of T[ri\) correspond to the Hugenholtz 
diagram , while the last two graphs of Fig. 4 (order f/^ terms of F [n] ) are equivalent to 
the last three graphs of Fig. 5 (order U"^ terms of /31^[Jo]) and correspond to the Hugenholtz 
diagram (^^^ reference 29). Therefore, we see that all the order f/^ terms cancel 

each other out in the effective action T[n]. 

If one were to expand the effective action (80) in powers of U, our approach reduces 
to performing the inversion method using as the expansion parameter. In this case, 
M//>i[Jo] in Eqs. (47-73) contains exactly all the order diagrams in W^[Jo] and can be 
expressed as Hugenholtz diagrams of order as well. Since all the Hugenholtz diagrams 
give value zero, all the derivatives of Wi (all the order diagrams) with respect to the 
density vanish as well. This implies that all the Ji>i vanish and consequently in the effective 
action the sum of Hartree energy and the exchange-correlation energy equals zero at zero 
temperature when there is only one electron present. 

The perfect cancellation of the Hartree energy and exchange-correlation energy means 
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that one is the negative of the other. This means that Yl'^i^i — —^n-oUon. Employing 
Eq. (77) for the most general case, we obtain 

Since Jo = Jo + Uon, we find that Jq = J at the single electron limit at zero temperature. 
The fact that Jq = J means that the final effective potential Si— ii — J-n is nothing but the 
lowest eigenenergy of the following single-particle Hamiltonian 

/l(x) = h l'ion(x) - H+J 

zm 

less the expectation value of J, which is exactly what one expected. 

D. Screening of Coulomb Interaction 

In the physical limit, n — i(f. In our formulation, lumps of charge fluctuation around the 
configuration i(fi interact with each other via 

UoV^K U =(U - UoDqo U) . 

As will be described in the discussion section, Uo(i(f)) — ib plays the role of the photon field 
here. Therefore, 

-00 (U - UoDoo U) 00 = (i0o U)oV-h(Uoi(l)) = (ib)o {U-^ - Do) o{ib) . 

Since U{x,y) = 6{tx — Ty)U{y.,y) and Qo^x^y) only depends on the relative time difference 
'Tx ~ Ty, we expect DQ{x,y) to depend only on r = r^; — Ty. Furthermore, since Qo{x,y) is 
antiperiodic in r^, Ty, and — Ty, one expects that ^o(^) y)Go{y, x) to be periodic in t^ — Ty. 

Introducing the spatial momenta p and q conjugate to the spatial variables x and y, we 
may write 

[U-' - Do] {u^, p, q) = j e'^^+'^^d^dy j^^ dT e^-^ [[/"^(x, y)5{T) - Qoix, y)go{y, x)] . 
When t/(x,y) = we find that U-^{x,y) = -^S{t^ - Ty)Vl6{x - y) and 



C27r)3 2,3 



where is the spatial volume of the system. Let us now write Do{i'n, p, q) as 
^o(^n,P,q) = j e'^^+'^^d^dy j^^ dT e""-^ go{x,y)go{y,x) . 
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In momentum space, oc q^5(p + q) and it is well known that this type of Coulomb 
interaction leads to infrared divergence (that occurs near q — > 0). In the presence oi —Dq, 
one may wish to calculate the zero momentum contribution oi —Dq. Let us define Do(m„) = 
Doi^n, P = 0, q = 0). Using Eq. (95), we find 

a,a' "^^ a 

Prom Eq. (101), we see that when x — x' , 
We therefore have 

where n — N^/ is the average electron density. 

Therefore at the stationary limit, where Vn — 0, the inverse propagator with nearly zero 
momentum behaves as 



q dn 



L 



3 



q^ + 47ie'^ — 

OfJ, 



analogous to the Thomas-Fermi results for electric charge screening. 

Except during the intermediate steps of computing the excitation spectrum, one deals 
with the time independent system. Note that Ane^dn/diJ, plays the role of in the screened 
potential e"'^'^ /r. This shows that the static interaction between charge fluctuations in the 
long wave length limit is a screened one instead of the bare Coulomb interaction. It should 
be noted that the use of this screened propagator dates back half a century. DuBois^^ re- 
placed the bare Coulomb interaction by the screened interaction in his study of electron 
interactions. Hedin^^ also used it to replace the bare Coulomb interaction in the expansion 
of the Luttinger-Ward functional,^^ resulting in the so-called GW approximation. The dif- 
ferences among the three mentioned approaches should be described. In both DuBois's and 
Hedin's formalism, they use the full polarization of the interacting system. The difference 
between their formalism is in the electron propagators employed. In DuBois's approach, the 
free electron propagator is used, while in Hedin's method, similar to the Luttinger-Ward 
functional, the full electron propagator is employed. In our approach, it is the polariza- 
tion of the KS system that enters the calculation and the electron propagator entering the 
diagrammatic calculation is the noninteracting KS Green's function. 
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E. Homogeneous Electron Gas 

Being the foundation for the LDA of the DFT, the homogeneous electron gas (HEG) is 
also a simplified model for a metal or a plasma. The HEG system is an interacting electron 
gas placed in a uniformly distributed positive background chosen to ensure that the total 
system is neutral. Due to the translational symmetry of the HEG system, the one-particle 
propagator functions only depend on the coordinate difference between two variables (space- 
time points) instead of on both variables. Consequently, each propagator (Green's function) 
carries a definite momentum even if the Coulomb interaction among electrons is fully taken 
into consideration. 

Investigation of the HEG system dates back to 1930s by Wigner,^''' who coined the term 
"correlation energy" to represent the ground state energy per electron after subtracting 
away the average kinetic and exchange energy. (The Hartrec energy is cancelled exactly 
by the interaction between electrons and the positive background, and by the Coulomb 
energy among the uniform positive charges.) Apparently, after choosing the Rydberg (the 
negative of the Hydrogen atom ground state energy) as the energy unit, one may express 
either the correlation energy or the ground state energy of the HEG system in terms of the 
dimensionless quantity = ro/r^, where = V being the volume of the system 
considered, being the total number of electrons, and = representing the Bohr 
radius. In fact, efforts have been invested to express the correlation energy as a power scries 
in Ts (and Vg In as well) at both the high density and low density limits. Lacking real 
experimental results to compare with, however, it is hard to assess how much improvement 
each increasing order of (or l/r^) can bring. 

Since oc I/tb oc e^, a small expansion is most naturally performed by treating 
the Coulomb interaction as a perturbation. Equivalent to the method of computing the 
vacuum amplitude, ^^'^^'^^ the time-ordered perturbative expansion developed by Goldstone^'^ 
formalized the Brueckner theory^^ and made a direct connection to diagrammatic expansion 
in calculation of the ground state energy. The electron propagator in diagrams under this 
formalism is that of the noninter acting electrons. An alternative formalism to compute the 
ground state energy (or the grand potential) is to utilize the full propagator (i.e., the one 
with self-energy included) in diagrammatic calculation. Under this alternative formalism, 
only two-particle irreducible diagrams contribute. Indeed, the work of Luttinger and Ward^^ 
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(as well as that of Klein^^) was exactly along this hne. Since the self-energy is not known a 
priori, this type of energetic expression is a functional of the self-energy (or the full Green's 
function), which must be determined via a stationary condition. 

An equivalent of the Brucckncr-Goldstone formalism^^'^^ was used by Gell-Mann and 
Brueckner^^ to compute the correlation energy of the HEG system. Given the long-range 
nature of the Coulomb interaction, it is not surprising that Gell-Mann and Brueckner iden- 
tified the occurrence of divergence as early as in the second order of e^-based perturbative 
calculation. To circumvent this unphysical divergence, they summed an infinite subset of di- 
agrams to arrive at a contribution proportional to Inr^. To eliminate unphysical divergence, 
DuBois''^ replaced the bare Coulomb interaction by the screened one. Starting with calcu- 
lating the vacuum amplitude, he expressed the ground state energy in terms of integration 
of the full electronic polarization over the Coulomb coupling strength, using a version of the 
Feynman-Hellmann theorem. While the full polarization is used, the electron propagator 
entering DuBois's formalism is the free electron propagator instead of the full propagator. 
Unfortunately, DuBois made a mistake in extracting the higher order terms of the ground 
state energy at the high density (small r^) limit. Although this error was later corrected 
by Carr and Maradudin,^^ according to Hedin,^^ this high density expansion violated at 
moderate values Ferrell's condition, which is based on the simple fact that the second 
order perturbation contribution to the ground state energy is always negative. 

Except for using the screened Coulomb interaction to replace the bare one and working 
directly at zero temperature, Hedin's approach is largely similar to the finite-temperature 
formahsm of Luttinger and Ward^^ in the sense that the full electron propagator is used 
as the fundamental variable. As a consequence, the higher order diagrams contributing to 
the full polarization appear different in DuBois's work and in Hedin's formalism. Instead 
of solving his own self-consistent equations, however, Hedin approximated the full Green's 
function within his formalism by the non-interacting Green's function to compute the ener- 
getics of the HEG system. Nevertheless, it is still possible to maintain the exactness within 
Hedin's approach if two-particle reducible diagrams arc properly included. We will provide 
an example of such a two-particle reducible diagram that would not be included in Hedin's 
approximate calculation but should be incorporated for theoretical soundness. 

Another important issue in obtaining the ground state energy of a many-electron system 
has to do with whether the finite temperature formalism (setting V ^ oo first followed by 
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T — )■ 0) or the zero temperature formalism (setting T — > and then followed hy V ^ oo) 
is used. In terms of diagrammatic expansion, there exist diagrams (termed anomalous 
diagrams by Luttinger et al^^'^"^) that are present (giving finite contribution) within the finite 
temperature formalism but are absent (giving zero contribution) within the zero temperature 
formalism. To be specific, a diagram is anomalous if within it there exist two electron 
propagators linking two different times (say Ti, T2 and of course Ti T2) in the opposite 
orders: ^(x, n, y, T2)^(w, T2, z, ti). 

Zero temperature methods typically rely on the Gell-Mann and Low theorem,^^ which as- 
sures that adiabatic transformation of the noninteracting ground state by gradually switch- 
ing on the interaction leads to an cigcnstate of the interacting system. ^^'^^ The energetic 
difference of this cigcnstate and the noninteracting ground state is what the zero temper- 
ature formalism obtains. This approach therefore makes sense only if the adiabatically 
transformed state is the ground state of the interacting system, which occurs only if the 
noninteracting Fermi surface is identical to that of the interacting system. For the HEG 
system, the perturbed Fermi surface remains spherical (identical to that of the unperturbed 
one) since the Coulomb interaction respects spherical symmetry and the background pos- 
itive charge distribution also has spherical symmetry. Consequently, for the HEG system, 
the anomalous diagrams should end up giving no contribution. Indeed, Luttinger et a/.^^'^^ 
illustrated how the contributions from the anomalous diagrams in this case are cancelled by 
the chemical potential shift. Luttinger and Ward^^ also showed how one may avoid anoma- 
lous contributions from appearing by expressing the grand potential as a sum of all possible 
hnked diagrams (including two-particle reducible ones). The key step there is to subtract 
from each self-energy part a number, which is given by that self-energy part evaluated at 
the Fermi surface. As we will illustrate later, under our finite temperature formalism it is 
not necessary to implement such an elaborate subtraction scheme because each two-particle 
reducible diagram is automatically accompanied by another appropriate diagrammatic sub- 
traction. 

What sets our self-consistent equation (78) apart from that of Luttinger and Ward and of 
Hedin is the variable to be solved for. It is the KS potential, instead of the physical Green's 
function, that enters our exact, self-consistent equation. In general, one needs to first solve Jq 
self-consistently using eq. (78) prior to the evaluation of the grand potential (or the ground 
state energy). When limiting to a homogeneous system with constant electron density. 
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however, Jq = const, becomes the only possibihty. Whether or not such a choice satisfies the 
self-consistent equation (78) depends on whether a uniform electron density truly represent 
the lowest energy configuration. For the present purpose of considering a high density HEG, 
we assume a constant Jq to proceed. Because a constant Jq can be easily absorbed into the 
chemical potential, the resulting KS Green's function carries a corresponding energy that 
consists of kinetic energy only. That is, rather than being an approximation as in Hedin's 
case, the single-particle Green's function carrying only kinetic energy represents exactly 
the self-consistent KS Green's function under our formahsm. Consequently, for the HEG 
system, the grand potential (or the ground state energy) may be calculated using eq. (80). 
In the remaining part of this section, we will first show how our formalism naturally avoids 
divergence and how one may use it to obtain the ground state energy of the HEG as a 
series in (and Inr^). We will then illustrate with an explicit example how the anomalous 
contributions are cancelled within our formalism, followed by a brief description of diagrams 
that would be missed within Hedin's approximation^^ in computing the ground state energy 
of the HEG system. To provide an easier comparison with existing results, we restore in this 
section the electron spins that have been suppressed thus far to simplify the exposition. 

In our definition of the energy functional Ey, see eq. (4), the —jJ^Ng term is included but 
the interaction among background charges is not included. Since most zero temperature 
formalism does not include the —/iNe term and for the HEG system the interaction between 
background charges is also included, the ground state energy in the literature will correspond 
to 

lim -r[n] + (iNe + -nbg- U-ribg 
in our formalism. Eq. (104) can then be used to arrive at 



Eg — hm 

/9-*-oo 



1 1 

-r[n] + iiNe + -nbg- U-n^g 



1 

2' 



Jn-n H — n- U-n H — nh,„- U-nh,s + lim 

2 /3->oo 



a=l 

y](£a - ^o) + lim 

a=l 



1 1 00 

-Trln (v^Ku) + -J2Uri] 



1 1 °^ 



(114) 



where Jq = Jq + U-n is used and = n for the HEG system is also employed. Note 
that the state label a — (p, a) includes both momentum and spin. For the HEG system, 
^ion + U-n — and thus eq. (94) leads to £q, = p^/2m -|- Jq. Therefore, the ground state 
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energy of the HEG system may be written as 

|P|<PF 2 1 fi oo ■ 

= 2 E Is + ^ 5^-' (*o '» ^)+J: r.W . (115) 

P y \_ i=2 _ 

where pp indicates the Fermi momentum and the factor of 2 in the first term of the right 
hand side comes from noting that there are two spin states associated with each momentum. 
Furthermore, the KS Green's function in this case may be written as 

1 g-iWn(Ta;-T;,)g-ip-(x-y) 

go{x,a;y,a') = —2_^— ■ , (116) 

where Ep = p^/2m is the kinetic energy of an electron carrying momentum p. By absorbing 
the constant Jq into the chemical potential in (116), one sees that the KS Green's function 
in the HEG system is indeed the free electron propagator and at the zero temperature, the 
new chemical potential (the original one subtracted by Jo) simply becomes p^p/2m. 

The diagrammatic expression of our first correction term Ti[n\ — |Trln('Do ^°U) — 
^Trln [I — Dqo U] is shown in Fig. 1. It can be seen that, upon being divided by the inverse 
temperature /3, our Ti[n] contains the exchange energy tx and all the ring-like correlation 
energy e' discussed by Gell-Mann and Brueckner^^ via the relation ri[n]//9 = Ne{ex + e'), 
with Ng being the total number of electrons. Since both Dq and U are diagonal in momentum 
space for the HEG system, we may write 

;^riW = ^El^[l-^o(g,^n)C/(g)] = ^/ T^,J2^n[l - Do{q,Un)U{q)] . (117) 

By comparing this equation to a derived result (eq. (30.16) of reference 48), it is also evident 
that our ri[n\/{/3Ne) indeed gives €x + e' of reference 54 as /3 — >■ oo. In Fig. 1, except for 
the first diagram, all other diagrams when evaluated individually exhibit infrared divergence 
due to the piling up of 1/q^ propagators.^^ To obtain the leading contribution, one must 
pay particular attention to the small |q| region. 

The polarization Dq = 6n{x) / S Jo{y) is defined in eq. (71). When the electron spins are 
included, n{x) — —^^goix,cr;x,a). Because the Coulomb interaction does not fiip spins, 
the KS propagator is diagonal in the electron spins and thus 

Do{x, y) = E ^0(2;, cr; y, (7)^o(z/, cr; 2 go{x, y)^o(y, x) . (118) 
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We therefore have 
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~l~ ^p+q ^p. 
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(27r)3 (£p+q - £p)2 + ^2 

where the third hne of the above equation is obtained using the technique described in 
section HIE, the fourth hne is obtained by rewriting the numerator of the third hne as 
— np(l — np+q) +np+q(l — rip), the fifth hne is obtained by a change of variable (p + q — >■ — p) 
in the second half and assuming the spherical property of rip, and the sixth line comes from 
the fact that ripfip+q is symmetric with respect to (p + q) p while the quantity inside the 
square bracket is antisymmetric with respect to (p + q) ^ p. Since £p oc |pp is a monotonic 
increasing function of |p|, while both rip = \l (^e^i^^-i^) + 1) and 1 — rip are nonnegative, we 
see that I?o(q, ^n) is a negative definite quantity. 

It turns out that for the HEG system one can further simplify the expression of Do(?, J^) 
by introducing a new variable u via v = u\c\\/m. We then have (by choosing q as the z 
direction in p and using the fact that rip = ^dpl) = n(p)) 
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As /5 — )■ oo, n{p) = 9{pp — p) and therefore dn{p)/dp = —S{p — pp), leading to 



m /3^oo (27r) 



Pf + TryPp + u - — )ln- — 

2q 4 (PF-f) +M 
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-ti tan h u tan 
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u u 

The exact expression (121) allows one to extract the limits of g — )■ oo and g — )■ 0, both of 

which are important for determining the convergence properties of the energy expansion. 
We have 

Do{q » 1, ^) -^^1^ + O {{q' + An-)-') , (122) 

L>o(g « 1, -) y [R,{u) + i?i(u) g2 + R^[u) q'] + O (g^) , (123) 

where 

R^{u) ^ PF-utaxT^— , (124) 

^iH = ~TT^7Z#^-^ ' (125) 



12(p^ + m2 
240 (p| + u2 



= • (126) 



With these asymptotic behaviors, we see from eq. (117) that Ti[n\ is finite. Methods for 
extracting coefficients associated with the e^-based expansion for e' can be found in references 
54 and 55. The e' = lim^_j.oo '^i[n]/{Nf,f3) — part contains in general r"-°lnrs and ri"-°^ 
terms when the energy is expressed in Rydbergs and expanded in power of (or e^). 

We now turn our attention to r2[n], whose diagrammatic expression is given in Figure 3. 
As mentioned earlier, within our formalism, rj>2[n] correspond to diagrams containing only 
Vq propagator and the KS electron propagator. Since the Vq propagator retains a finite value 
even when its associated momentum approaches zero, each of the diagrams corresponding to 
rj>2['^] takes a finite value. The absence of divergence no longer holds when one perform high 
density ( small r^) expansion, as we will illustrate explicitly later. To make easy comparison 
with existing HEG studies, one employs Dyson's equation, 

= U + UoDqoVq (127) 

<ww = • ♦ -)- 
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to decompose the propagator T>q (of order and higher) into the sum of a bare Coulomb 
propagator U (of order e^) and UoDqoVo (of order and higher). The reason that one 
should not expand further and write "Do = f/ + f/ oDqo U + U oDqo U oDqoVq is because the 
UoDqoU term causes infrared divergence due to the momentum integral J dq/q^. In fact, 
this is exactly what causes (in the second diagram of Figure 1) the infrared divergence, 
motivating the summation of the ring diagrams. 

To illustrate the main points, let us begin by examining the first three diagrams of T2[n] 
and employ the decomposition rule mentioned above: 




Let us start with the diagrams to the left of the equal signs in (128-130). We will for now 
keep the bosonic propagators general, that is, allowing them to carry frequencies. After 
that, we will discuss the (a) diagrams in (128-130), followed by the (b) diagrams and then 
the (c) diagrams. 

In (128), let the vertical boson propagator, denoted by Bi, carry momentum q (upward) 
and frequency u. Let the horizontal boson propagator, denoted by B2, carry momentum 
q' = — (pi + P2 + q) (leftward) and frequency u'. Using techniques described in section HIE 
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for frequency summation, one obtains the following generic result 
/3V 2 ^ f dpi dp2 dci -Bi(q,i/) S2(q',i/') / Up^ - Up^ 



E 



^ P"^ , J (27r)3 (27r)3 (27r)3 -ii/ + £p^+q - £p^ ii/ + £p2+q - |^£p2 -£pi + 
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^^Pl+q' - £^pi - il^' £p2+q' - ^P2 + ^P2+q " ^Pl+q + ^(^ " J 

or equivalently, 

PV 2 ^ f dpi dp2 dq Bi{q,u) B2{q',u') ( rip^ - Up, 
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4 /S^ 7 (27r)3 (27r)3 (27r)3 -ii/' + Sp^+q' - Sp^ iu' + Sp^+q' - £p2 I £p2 - ^pi - ^{^ + 



Tbp^ JT-pi+q ^ Tbp^ ^P2+q _| ^P2+q ^pi+q I (132) 



Spi+q Spi ^P2+q ^P2 ^P2+q ^pi+q ~^ ^(^ ^ ) 

Note that the factor 2 associated with ^ comes from the two possible spin states for an 
electron. 

In (129), let the boson propagator on top, denoted by carry momentum q (towards 
lower-right) and frequency v. Let the boson propagator at the bottom, denoted by i?2, 
carry momentum q' = (pi — P2) (towards upper- left) and frequency v' . Using techniques 
described in section III E for frequency summation, one obtains the following generic result 

^^v^ [ dpi dp2 dq Bi{q,u) B2{q',u') 

2 fi^^J (27r)3 (27r)3 {2^f -^v + £p,+q - £p, -^v' + ep^^^ - Sp,,^^ ^ ^''^^^''^ ''^^^''^ 

^pi+q ~ ''^P2+q _|_ ^pi ~ ^pi+q | ^pi ~ '^p2+q I (133) 



—iu' + £pi+q — £p2+q + ^pi+q ~ ^pi ^(^ ^ ^0 + ^pi ~ ^P2+q 

In (130), let the boson propagator at the right hand side, denoted by Bi, carry momentum 
q (downwards) and frequency u. Let the boson propagator at the left hand side, denoted 
by B2, carry momentum q' (also downwards) and frequency u'. Using techniques described 
in section III E for frequency summation, one obtains the following generic result 
l3V 2'^ -1 ^ f dpi dp2 dq dq' Bi{q,i^) B2{q!,u' 



2 2/3 / j§^np{l -np)j^J (27r)3 {2ny (27r)3 (27r)3 -zu + £p,+q - Cp, -lu' + £p,+q, - Sp, 

\ /3^^pi(l - ^pl)^P2(l - ^P2) + ^^Pl(l - ^Pl)— 7 



^P2+q' ^P2 
+ ^P2+q' ^P2 

+/3np,(l - Up,) ~ ""P^ + ^P^+q' - ""P^+q ~ 1 , (134) 

-IV + ep^+q - Ep^ -iiy' + £p2+q/ - ep2 -lu + Cp^+q - Sp^ J 

where the factor —1/ |2/3 / (^^^p(l ~ ^p)| — B>Q^{q — 0,1/ — 0) arises from the fact that 
in a translationally invariant system, momentum conservation at each vertex demands that 
the inverse density correlator must carry zero momentum and frequency. 
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For the diagram (a) on the right hand side of eq. (128), -Bi(q, z/) = 47re^/q^ and 
i?2(q',z^') = 47re^/(pi + P2 + q)^ are both frequency independent. Thus, one may sum 
over both v and v' . And upon doing so, we obtain 
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If we divide the quantity above by A^e and then take the zero temperature hmit, it gives rise 
to (using as the unit for momentum) 
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Rydberg, 



which is exactly the e^^"* term of Gell-Mann and Brueckner.^^ 

Diagram (a) on the right hand side of eq. (129) is one of the anomalous diagrams^^'^"^ that 
give rise to finite contribution as the T — )■ hmit is taken within finite temperature formahsm 
but are absent within zero temperature formahsm. For this diagram, i?i(q, i/) = 47re^/q^ 
and -B2(q', = 47re^/(pi — pa)^ are both frequency independent. Thus, one may evaluate 
its zero temperature contribution by summing over both u and u'. Upon doing so, one 
obtains from diagram (a) of (129) the following anomalous contribution 
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dpi dq dq' 47re^ 47re^ 

7T^P^Pl(l — '"'Pi)'"'Pi+q''^pi+q ) 



(27r) 



q" iq 



(135) 



where the last expression is obtained via the following change of variables: pi(2) — Pi(2) — q, 
eliminating p2 by q' = pi — P2, and q' — )■ — q'. This expression can be further simplified via 
the following definition 



leading to 




/(P) 



-\/(47re 
-1/(47re^ 



dq 1 



(27r)3q2''P+'i' 
dpi 



{2n) 



/3np,(l - ripj/ (pi 



2 2 f dpi dup^ 



(136) 



(27r)3 9/i 

As mentioned earlier, each two-particle reducible diagram within our formalism is ac- 
companied by a corresponding diagram that will eliminate anomalous contribution when 
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applicable. Diagrams in eq. (130), appearing only under the effective action formalism, are 
such diagrams. They are neither present in the zero temperature formalism of Goldstone^° 
and Brueckner^^ nor in the finite temperature formalism of Luttinger and Ward.^^ Dia- 
gram (a) on the right hand side of (130) corresponds to the case i?i(q, z/) = 4'7re^/q^ and 
i?2(q', i^') = 47re^/(q')^, both being frequency independent. Thus, upon summing over both 
u and u' one obtains 




2{-lfV{47rey 

dp 



dpi dp2 dq dc( Aire^ Aire^ 



(2vr) 



12 



2/3/(|g.^p(l-^P 
y(47re2)2 [ f dp dup^^^^ " 



g2 q'^ 



r dp drip 
J (27r)3 ^^l 



(27r)3 9/i 



(137) 
(138) 



When combined with eq. (136), one obtains 





-V{ATxe 



2\2 



dp dup 

(27r)3 dfi 



{f-fY, (139) 



where the overline symbol is defined as / = / (^^^)f/(p)- Apparently, the expression 
(139) is in general negative unless / = /• We will show that this happens only at zero 
temperature and only if the system has spherical symmetry. 

In order to have /(p) = /, p can only have support at a constant /(p) surface. This 
is achieved when T — )■ 0, where Up — )■ 9{fi — Sp) and dup/dfi = 6{fi — Sp) forcing p to lie 
on a constant Ep surface. Furthermore, /(p) = / demands that /(p) depends only on the 
magnitude of p, i.e., /(p) = f{p). This can only be achieved if Up depends only on |p| = p. 
For the HEG system, = p^/2m, thus Up = n{p) and fi = p|,/2m. We therefore have 
dup/dfi = m5{p — Pf) Ipf^ fixing the length of p. Now rip+q = 0{fi — (p + q)^/2m)||p|=p^ = 
0(— cos-i? — 2^)) with being the angle between q and p. Thus, in the integral defining 
/(p), although p defines the z direction of vector q the integral is independent of the choice 
of p. Therefore, the anomalous contribution (136) may be written as (when T — t- 0) 



dpi m 

W¥p^ 



S{Pi-Pf) 



dx 

1 ^0 



q dq Aire 
(27r)2 g2 



X 



2pi 



(140) 
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while the corresponding subtraction term (138) may be written as 



V 



I 



dpi m 



V 



(ipi m 

Kpi - Pf) 



Kpi-Pf) 



1 roo ^2 

dx 
1 Jo 



(27r)2 g2 



dx 



g2 dq 47re^ 
(2^^ 



ei-x-^) 



:i4i) 



cancelhng exactly the anomalous contribution (140) in the HEG case. 

Within the framework of Luttinger et a/.,^^'^^ the anomalous contribution is cancelled 
by the chemical potential shift, the difference between /x(T — )■ 0), under finite temperature 
formalism, and /i(T = 0) = p\/2m, under zero temperature formalism. Within the current 
finite temperature framework, however, the /i(T — )■ 0) is identical to /i(T = 0) and the 
cancellation of anomalous contribution is explicit. As we will show later, however, diagram 
(b) of (129) can't be cancelled by diagram (b) of (130). This is because diagram (b) of (129) 
is not an anomalous diagram. 

Before moving onto (b) diagrams in (128-130), let us remark that diagrams within our 
formalism, i.e., diagrams on the left hand side of (128-130) with -Bi(2) T^o yield only finite 
contribution. The decomposition made in (128-130), however, may introduce divergence as 
we will illustrate using the (b) diagrams that correspond to the main results of DuBois.^^ 

For the (b) diagram of (128), i?2(q',«^') = 47re^/(pi + P2 + q)^ while 5i(q, z/) = 
(47re2)2L)o(q, z^)/ {q^[q^ - 47re2Do(q, v)]}. Therefore, one may sum over z/' to simplify the 
expression. Using again the methods described in section HIE for frequency summation, 
one obtains 




V (47re' 



2\3 



4 (27r)9 13 



V / rfpirfparfq— — 



^o(q, ^) 



47re2Do(q, i^)) (Pi + P2 + q)^ 



n 



pi+q '^pi 



^P2+q ^P2 



K(W)3 

^ (27r)9 



-If + £pi+q £pi If + £p2+q £p2 

^o(q, 



V / dpirfp2c/q^— 



2-47re2Do(q, i^)) (pi - P2) 



X- 



^pi+q '^Pi 



n 



P2+q 



n 



P2 



_. , _ _■ , _ ' (142) 

tU + £pi+q £pj + £p2+q £p2 

where the last expression is obtained by changing variables: P2 + q — >■ — P2 followed by 
P2 ~^ P2- For the high density expansion, where is treated as a small parameter, the 
major contribution in (142) comes from the region q — )■ 0. One thus writes 



n 



p+q 



n 



p I ' l^p+g 

IqKO 



-/3np(l - Tip)(ep+q - ^p) -q cos 1) 6 {p - Pf) 
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where ■& is the angle between p and q. As T — )■ 0, | F{iyn) — )■ f^-^(^) if -^('^) does 
not have pole strength greater than one. Making a change of variable u = u\q\/m and 
treating g as a small quantity, the major contribution of the (b) diagram of (128) is given 
by 



-co 

"I, 



du f m6{pi — pf) cos^i m6{p2 — Pf) cos^2 

— / dpi dp2 n a ^ 

2tt J —iu + pFCOSVi iu + pF cos V2 

Anqdq- ^^^^'^^Z^) 



Pi-P2)^yo q"^ - ATce'^Do{q,uq/m) 

(Ane^Y /"°° '^'^ cos -i?! c/ cos "i^i cos 'i?2 cos ■(92 



\^ ffl — -— / / X 

(27r)9 J _^ 2-K J _i —iu + pf cos i^i —iu + pFCos'd2 

/'4.,d„ ^. (143) 



I COS -i?! — COS "(92 1 7 q^ ~ ATie'^DQ{q,uq/m) 

Apparently, the 1/| cosi^i — cos'(92| factor in the integrand causes an undesirable divergence 
in the d integral that is absent if one does not decompose the diagrams using (127). A finite 
result emerges only when one combines (143) with the (b) diagram of (129), which we soon 
turn to. 

Before elaborating on the evaluation of the (b) diagram of (129), we wish to point out that 
this is the type of diagram (two-particle reducible ones) that was missed in Hedin's approx- 
imation'^^ but should have been included for the exactness of the theory. For this diagram, 
52(q',z/') = 47reV(pi - P2)2 while i?i(q,z/) = {A7ie^fD,{^,v)/{ci^[^^-Ane^Do{c^,v)]}. 
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Summing over v' via methods described in section III E, one obtains 




= V 



"^^^ "I" £^pi ^pi+q 



X 



p (27r)9 (pi-P2)^glg' 

[nnnA-n — ''^m ) ( ''^ni +n — ''^n, ) 2p?2ri, ( 



^P2+q(^Pl+q ^Pl;* _ 
"I" ^pi+q ^pi, 

1 Doiq,i^) 

'2 - Ane'^Do{q, u)) 

l^P2+q ~ ''^P2)(^pi+q ~ ''^pi) _j_ 2/3?2p^ (1 — "n-pjrap^ 
+ £pj+q -Sp^Y -iv + epj+q - e 



^ i-Vm) 



'TT 



'Pi. 

/•°° COS 



— tt^ -r tpi+q — tpi J 
l?! COS ^2 d COS i9i d COS '(92 



ist/i — cos'iy2| J q'^ — A-Ke^DQ^qjUq/r 

(47re2)^ /"°° '^^ /"^ '^t^^Pf dcos-d f dq DQ{q,uq/r.^^ 
(o.rr\3 I 27r J _i -iu + Pf COS 'd J {2iTy q'^ ~ Ane'^Do{q,uq/m) 

(144) 



X 



/m) 



' —00 



(2vr, 
■ q''^dq'dcosd' 
IJ (27r)2 (q')2 
= (-Fm)Li + (2\/)L2 , 



/m) 



cos'd' — 



' ) 



2PF 



where the first term inside the square brackets after the second equal sign is obtained by 
adding to it an equivalent expression with Pi + q — !■ — Pi, P2 + q — ^ —V2, ^ and 
then taking the average, while the second term inside the same square brackets results from 
changing u — t- —u. The L2 part, where the dummy variable of the integral is switched from 
P2 to q' = P2 — Pi) can be cancelled by one of the terms contributing to the (b) diagram of 
(130). The Li part, however, upon taking the T — i- and g — )■ limits, may be combined 
with (143) to yield a finite expression 

Vmpp (4716^)^ du cos T?! d cos cos -(^^i cos -(99 



du 

2 {2ny~J_^2^ 
cos -di — cos {^2 



271 d{q 



iu + pf cosdiY — m + Pi? cos'i92 
2^ Do{q,uq/m) 



(145) 



costal — cos i?2 1 y q'^ — 4:ne'^Do{q,uq/m) 

Note that this expression, in agreement with Carr and Maradudin,^^ carries a different sign 
when compared to the original result of DuBois.^^ 

It was conjectured before^*^ that cancellation of diagrams of the following form always 
hold true for HEG 



O * ^3 
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Here, black circles denote parts of the diagram that are connected to each other via two 
propagators. If this is the case, then the contributions of the (b) diagrams of (129) and (130) 
will cancel each other. We don't expect this to happen since the (b) diagram of (129) is not 
anomalous. To illustrate that the (b) diagram of (130) does not eliminate the (b) diagram 
of (129), we now proceed to evaluate the (b) diagram of (130). 

Here, 52(q', z/') = 47reV(q')^ while Si(q, u) = {AneyDoici, z/)/ {q2[q2 - Ane^Doiq, u)]}. 
Summing over u' via methods described in section HI E, one obtains 
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J2 J dpidp2 rfqrfq'^^^ 
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(2tt)^ m q^i^q^ — Aire'^Do^q, uq/m)) 
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'P2y 
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(q')2 



6{fi — epjm/q 5(/i — epjm^/g^ qpp cos-di 



-iu + Pf cos-i?! 
-(2\/)(47re^)V(27r)3 

2/(|^%-^p) 



-2M + Pf cos-i^i) 
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Do{q,uq/m) 
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(2vr) 
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(27r)3 g2(^g2 _ /^Txe^DQ^q, uqj m)) 
rfq' 0(-cos^'-^) 



(27r) 



(qO 



5(Pl 



PF)Tn/pF ^ 5{pi — pir)mcos'i9i 



— m + Pi? costal {—iu + Pf cos-diY 



i2V) 



2\3 



f27r) 



du Art mp F d cos § I 
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-iu + Pf cos di 
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Pi? COS'(?i 



—iu + Pf cos-i^i 



X 



dq 



Do{q,uq/m) 
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2pF ' 



(27r) 



iq 



(27r)2 (f — 47re^Do('?; uq/m 
The contribution of (146) can be easily divided in two by explicitly expanding the two parts 
inside the round parentheses. It is obvious that the contribution associated with 1 cancels 
exactly the {2V)Lo part of (144) while the contribution associated with pf cos i^i cannot 
cancel the {—Vm)Li part of (144). 

For the (c) diagrams of (128-130), both Bi and B2 are of order (47re^)^, leading to 
contributions of order {e^Y ^^^l higher. The (c) diagrams thus already lead us beyond what 
was studied by DuBois"^^ and by Carr and Maradudin.^^ The last two diagrams of r2[n] 



(146) 



62 



(shown in Figure 3) gives rise to the E'^ term of Carr and Maradudin^^ when the Vq hnes 
are each replaced by U, the first term in the decomposition of Vq. In principle, one may 
go on to study terms of order (e^)^; we will not, however, delve into this endeavor since 
this is not the primary aim here. We would like to emphasize the following points. First, 
unlike conventional based perturbation theory, the formalism presented here naturally 
avoids divergence. This is shown by the fact that each diagram in our formalism contains 
no singularity while attempts to perform based purturbation necessarily require further 
re-grouping, such as combining the (b) diagrams of (128) and (129), to tame the divergence. 
Second, even if one were to pursue expansion, using Dyson's equation (127) within our 
formalism still makes the task straightforward. In fact, as shown in this section, the To[n] + 
ri[n]+r2[n] part when applied to the HEG already contain the celebrated results of Carr and 
Maradudin.^^ Third, the removal of anomalous contributions that required special attention 
within the formalism of Luttinger et aL^^'^^ becomes automatic under this formalism. 

V. EXCITATIONS 

To obtain information regarding the excitations, one needs a time-dependent probe, al- 
though one with infinitesimal amplitude is sufficient. Runge and Gross^^ extended the 
correspondence between the external probe potential and the ground state charge density of 
the DFT to the time-dependent case. This relationship provides the foundation for studying 
excitation energies under DFT. Fukuda et al}^ expressed the excitation energy condition 
using effective action formalism, without explicit connection to the Kohn-Sham formalism. 
By introducing time-dependent Kohn-Sham orbitals while assuming time idependence of 
the orbital occupation numbers, Casida^*^ derived via linear response theory a self-consistent 
condition on the density matrix response that leads to determination of excitation energies. 
Along a similar line, Petersilka et al.^^ proposed the so-called optimized effective potential 
(expanded in time-dependent Kohn-Sham orbitals) to tackle the problem of excitation ener- 
gies. Since there exist formalisms to extract the excitation energies of the system provided 
that the UDF is known, our result for excitation energies should not be considered novel. 
The reasons for this section are twofold. First, we would like to exphcitly show that the 
excitation energies can be obtained using the formalism of section III without introducing 
time-dependent orbitals. Second, although it is possible to find derivations of formulas^° 
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similar to those that will be shown here, they do appear slightly different and thus a self- 
contained exposition may be helpful. 

Intuitively speaking, by varying the frequency of the probe, one seeks the fre- 
quency/energy where the amplitude of the response function diverges. Indeed, it is known 
that the spectral representation of the correlation function has poles at the excitation ener- 
gies of the system^^. Having obtained the effective action r[n], we also note that the second 
(functional) derivative of r[n\ with respect to the local electron density is the inverse of 
the density-density correlation function. Therefore, any pole associated with the correlation 
function becomes a root of the effective action. It can be shown that upon analytic con- 
tinuation the correlation function, obtained using the imaginary time (finite temperature) 
formalism, can be turned into the response function of the real time. Below we briefly il- 
lustrate this point. Readers interested in more details can find an extensive exposition in 
reference 48. 

Prom Eqs. (13) and (28), we know that 

S{PW[J]) 



n[x] 



and J{x) 



5J{x) 

sr[n] 



Sn{x) 

One then considers 

_ , SJ{x) _ f 6^T[n\ 5n{z) _ f S^Tjn] S\f3W[J]) 

5J{y)~ J 6n{x)5n{z)SJ{y)~ J Sn{x)5n{z) 6J{z)5J{y) ' ^ ' 

Although a time-dependence of J is introduced to probe the excitations, in the end we will 

return to a time-independent source {J{x) — )■ J(x)) while computing the excitation energies. 

As we will show below, it is most convenient to go to the zero temperature limit to compute 

the excitation energies. 

Note that 

Ij{x)5^J^) ^ ~ [^(^)^(^)]t + [^(^)]t Hy)]T = - Hx)h{y)]T + n{x)n{y) 

= - [{n{x) - n{x)){n{y) - n{y))]r^ = - [n{x) n{y)]r^ (148) 

where n{x) = ■ip''{x)ilj{x) is the electron density operator and the square bracket []t indicates 
an (imaginary) time-ordered thermal average such that 



Tr 



e-WlT {^^{h)62{t2))] Tr [e-^T (01(^1)02(^2))' 



It Tr[e-^H[j]^ Z[J] 
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Let us denote by {\i)ex}'^o the eigenstates of the Hamiltonian H[J] with the corre- 
sponding eigenenergies {£e}'^Q- The spectral representation is obtained by first writing the 
time-ordered product (assuming that operators Oi(ti) and 02(^2) are bosonic) as 

e.(^|r (di(h)d2(t2)) \e)ex = ^(^1-^2) ex{£\di(h)d2(t2)\£)ex+0{t2-h) ex{Wt2)di(h)\£) , 

then by inserting the identity operator \i')ex ex{£'\ between the two operators, and finally 
by multiplying by e*'*'^*^'"*^) and then integrating over the time variable ti — t2- Proceeding 
in this way, one will then obtain information on S^/ — £1. Since knowing Q^/ = S^i — £q also 

provides complete information on E^i — 81, one may also focus on £ = by taking the limit 
(5 ^ 00. Let 

and 

-iy(^)(x,y,i^„)^ rd(r.-r,)e^-(--^)iyP)(x,y)^ T ^t^^ , (150) 

Jo J -00 ^^n - 

with Vn — 2Tm/ P and 

A{u) = e^^^'^Y.^-^'' (^"''" - 1) + - e.(^|^(x)|m)e.e.(m|n(y)|£)e. . (151) 

Note that n(x) measures the deviation from the thermally averaged electronic density nr(x). 
Since the expression (149) is evaluated at static J(x), the Hamiltonian contains no time 
dependence. We may thus write n{x) — e^'^'^n(x)e~^'^"^. 

Let's now express using real time the retarded correlation function (7^), also called 
the response function, and the advanced correlation function {A) as follows (with n{x) = 



A{x,y) I \ ie{ty-Q 

-i9{tx - ty) 

i9{ty - tx) 



J2 e^mJ]-e^) [e^(*-*^)(^^-^-),,(£|n(x)|m)e.e.(m|n(y)|£)e. 
_e-i(t.-t.)(£,-£^)^^^^|^(y) |m)e. e.(m|n(x) |£)e.] • (152) 
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Taking the Fourier transform of the response function, we consider 



J2 e-^'^ ex(^|n(x)|m)e. e.{m\h{y)\i), 



doj' A{uj') 



oo 



(153) 



27r uj — uj' ± ir] 

Comparing Eqs. (150) and (153), one finds that substituting ii/n oo + ir) {ou — irf) in the 
imaginary time-ordered correlation function leads to the retarded (advanced) correlation 
function. The validity of this analytic continuation was discussed by Baym and Mermin^^. 

As T ^ 0, e-^^[-^l = e-/^^o [l + C>(e-'^(^i-^«))] , and e^^M = e^^^ [l - 0{e-^^^^-^°% 
Under this limit, we may rewrite the spectral weight A^uj) as 

hm = V(e-'^(^--^°)-e-'^(^^-^°V('^ + ^^-^-)-(^I^WIWe.e.(m|n(y)|£)e. 

t,m 

^0(g-/3(£i-£o)) 
= [S{u + Ee-Eo)eMf^{y)\^)exeMn{^Mex 

e. 

-5{u - {Et - £:o))e.(0|n(x)|£)e.e.(^|ri(y)|0)e.] 
= + E,- £o)n,%(y)n,,e(x) - 5{uj - {E, - £o))n,* e(x)n,,,(y)] (154) 

where n^,e(y) = ex{(\n{y)\^) ex ■ 

When continued to the retarded correlation function (response function), the density 
correlation function VT^^) reads 

'n:,(x)n£,e(y) n;g(y)nf,e(x) 



lim I^(^)(x,y,c^) = J] 

I 



(155) 



CO — fl£ + if] CO + fli + if] 
with 

ne = Ee-Eo. (156) 

Provided that the amplitudes exi^l^li) ex are nonzero, we see from Eq. (155) that uj + ir] — 
± {El — Eq) are simple poles of H^*^^)(x, y,a;). Furthermore, we also see that 

(W^(2)(x,y,a;))* = W^(2)(x,y,-a;*), (157) 

and when ui is real 

(M^(^)(x,y,a;))* = iy(^)(x,y,-a;). 
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Let us also define 

r(^)(x,y)^ 

' 6n{x)6n{y) 

Eq. (147) may thus be rewritten as 

-S{x -y) = J dz r(^) {x, z) W^"^^ {z, y) . (158) 

Since eventually, must agree with Ty in the equation above, if > r^, we must have 
> Tj; as well. Similarly to Eq. (150), the Fourier transform for F^^^can be written as 

F(2)(x, z, iun) = / d{T, - T,) e^-"(^^-^^) F(2)(x, z) . (159) 
./o 

The inverse transform of (150) and (159) can be written as 

n 

T^^\x,z) = i^e-''^"(^^-^'")F(2)(x,z,^i/„) 



and 

f/3 



^rfr^y^z F(2)(x, ^)iy(2)(z, y)^ jdz^Y. e-^'^"(^^-^^)F(2)(x, z, -ii.„) iy(2)(z, y, ii.„) . 

Since 5{tx - ^v) = ;| Z^n e"'''"^'^'^"'^^'^ one obtains 

-5(x - y) = y"dz F(2) (x, z, -ii/„) (z, y , ii/„) = ^dz M^^^) (x, z, -ii/„)F(2) (z, y , ii/„) , 

(160) 

which under analytic continuation becomes 

-5(x-y) = Jdzr^'^\x,z,-uj)W^'^\z,y,uj) ^ Jdz W^'^\x,z, -uj)r^'^\z,y,uj) . (161) 

From Eqs. (157) and (161), one has 

(r(2)(x,z,a;))* = r(2)(x,z,-a;*). (162) 

Multiplying the LHS and the middle expression of Eq. (161) by {u ^ Qe + irj) and then 
setting CO ^ ±fli — irj, we see that 

Jdzr^^\x,z,-uj ^ -ne-iri)nl^{z) = 0, (163) 

and y"dzF(2)(x,z, -uj +ne - irj) n£,e(z) = . , (164) 
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That is, (y) become eigenvectors of F^^^ (x, y , —u) with zero eigenvalues. The matter of 
finding excitation energies and the corresponding electronic densities thus reduces to finding 
for r^^) the eigenvectors with zero eigenvalue. Since we are interested in obtaining the ex- 
citation energy under the physical condition, J(x) — ?> 0, this also means that the derivatives 
of r above are evaluated at the ground state electronic density in the zero temperature limit. 

Prom the expressions (65) and (70), one sees that the effective action is spht into the free 
particle part Fq, the Hartree functional and the exchange-correlational functional 

r[n] = ro[n] + ^no Uon + T^eN = ToN + r,^t[n], (165) 

where F^cW — Yli'^i^ib^]- Letting A(x) be the eigenvector, the excitation condition be- 
comes 

y"dzr(2)(x,z)A(z) = 0. (166) 

After splitting the effective action into Tq and Tint the eigenvalue equations (163-164) may 
be expressed as 



- JdzT^^\y,z,-uj)A{z) = JdzT[^{y,z,-uj)A{z) . (167) 

Multiplying both sides of (167) by WQ'^\x,y,(jj) and then integrating over dy, one obtains 

A(x) = J dydzWi'\^,y,co)TS:l{y,z,-u)A{z) . (168) 

Note that Wo[Jo] describes our Kohn-Sham system, a constructed non-interacting system 
that produces the same ground state electron density as that of the physical system consid- 
ered. From Eq. (155), one can write down W^g^^(x, y, w) in terms of the excitation energies 
associated with the Kohn-Sham non-interacting system: 

e 

where 



uj — uji + ir) uj + uje + irj 



(169) 



uji = Ee-Eo, (170) 

and E£ is the energy of \i)ks: the £th state of the many-particle Kohn-Sham system with 
Jo(x) chosen to generate the correct ground state electron density. Therefore, for any £, 
is simply the sum of single-particle energies Note that in Eq. (169), n£(x) is defined by 
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ks{^\n{:x.)\0) ks except that \£)ks now describes the £th state of the Kohn-Sham system, not 
the physical system considered. 

We seek a general solution for A(x) of the form 



A(x) = ^ [aene(x) + 6^ n|(x)] 



(171) 



Evidently, if a frequency ou leads to a solution {6^}, then —a)* should lead to a solution {a^}, 
which plays the role of {b}}. Substituting Eqs. (169) and (171) into Eq. (168), we find that 

^ [aene{x.) + ben}{x)] 

n^(x)n^(y) n^(x)n^(y) 



X] \dydi 



u — uj£ + ir) u + uj£ + if] 



r2(y, z, -oo) \avnv{^) + 6,.n^(z)] . (172) 



Equating the coefficients associated with n^(x) and n|(x), we find that 



X f dydz— ^^^^i^i-7--r2t(y,z, -uj) [a^/n^/(z) + bem},{z) 
(1 J ' 



u — uJi-\- irj 



Let us define 



Ke,e'{uj) = J dydznl{y)rlll{y,z,-uj)nl,{z) , 
and we obtain the following matrix equation 



(173) 
(174) 



= (a; + irj) 



-1 
1 



A 
B 



;i75) 



Y(a;) K(a;) A 
K*{-u*) Y*{-uj*) j 

where (^4)^ = and {B)^ — b^. Evidently, one seeks uj such that 

Y(c:;) K(c:;) 
K*{-u;*) Y*{-uj*) 

As mentioned earlier, one anticipates {A(u;)) — (S*(— a)*)). To see this, we perform the 
change ou —ou* in (175) and find that one can then rearrange the resulting equation into 



det 



(a) + irj) 



-1 
1 







(176) 



Y{uj) K{uj) 



B*{-uj*) 



K*{-uj*) Y*{-uj*) / \ A*{-uj*) 



{id + 17]) 



-1 W B*{-u*) 
1 j \ A*{-u*) 



;i77) 



69 



which is identical to Eq. (175) except with a)*)) playing the role of {A{u;)) and 

(A*{-uj*)) playing the role of {B{uj)) . 

We now compare Eq. (175) with similar existing results. In references 63 and 20, equa- 
tions similar to (175) were obtained, and those will be identical to Eq. (175) provided that 
K*(— a;*) = K*(a;). This will happen if r(x, y, — w) = r(x, y,a;) for real cu. 



Below, wc will obtain the effective action using a classical variable itpc that corresponds 
to the saddle-point of the auxiliary field path integral. At the physical condition J = 0, 
i(fc is interpreted as the electron density of a self-consistent Hartree solution. Our Hartee 
problem is not of the conventional type, but rather similar to what Kohn described in his 
Nobel lecture.^ In the conventional Hartree calculation, the wave functions obtained may not 
be orthogonal to each other due to the fact that each particle's wave function is solved with 
a different potential. In the method below and mentioned by Kohn^, the electric potential 
experienced by every electron is the same. Another difference between the method below 
and the aforementioned Hartree methods^ '^^ is that the integral of the Hartree density in 
our method is not necessarily an integer due to the possibility that the density correction 
term may have a nonzero integral. 

The saddle-point method below is quite different from what was described in the pre- 
vious sections. First, although the diagrams in the saddle-point method are all connected 
diagrams, they are not one-particle irreducible (IPI). Second, unlike the method presented in 
previous sections, the computation of the effective action now requires no further functional 
derivatives of /3H^[J] with respect to J evaluated at Jc, while in the formalism mentioned 
in previous sections, one needs to compute higher order derivatives of ^VTfJo] with respect 
to Jo (see Eqs. (67-73) ). 



VI. SADDLE-POINT AS AN ALTERNATIVE FORMALISM 



A. Evaluation of e-^^^^t-^l via expansion around the saddle-point 



The path integral (23) 
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may be evaluated by the saddle point method. The extrema condition gives 

51[<t>\ 



5(t){x) 



— —iJ{x) . 



Since 



SI[cl>] 



5(l)(x) 



{Uo^,)^ - Tr 



6(l){x) 



(j>—>-ipc 



{U°^c)a,- J dy gc{y,y) {iU{y,x)) , 



(178) 



where Qc{y,y) = G^^^^{y,y), we obtain with U{x,y) = U{y,x) 

J{x) = {Uo{iip,))^ + Jdy U{x, y)g,{y, y) ^ j dy U{x, y) [iip,{y) + g,{y, y)\ . (179) 

As J — > (the physical condition), when U is invertible such as Coulomb interaction, one 
must have i(pc{x) — Scix, x). Note that the negative of the diagonal element of the Green's 
function —Gdx, x) is the particle density corresponding to the following Hamiltonian 

V2 



dx i/j^x) 



2m 



+ ^^ion(x) - IjL + i(Uo(pc)x 



V'(x) . 



The saddle-point equation therefore produces a Hartree-like equation: the Green's function 
depends on the input particle density i(pc{x) and is required to produce the same particle 
density i(fic{x) in the end. 

When J 7^ 0, one can still view (179) as a generahzed Hartree equation in the following 

sense. Remember that the inverse Green's function Qc^i^^y) given by (with 6{x — x') — 
5(t-t')5(x-x')) 



and may be rewritten as 



^T- ^ ^ion(x) 

2m 



b{x — x') , 



(180) 



dr- h 'Uion(x) + J{x) - A* + UoUh 

2m 



5{x - x') . 



(181) 



That is, we now view the potential as given by Vion(x) + J{x) and the Hartree particle 
density nuix) — i^c{x) — {U~^oJ)x. This interpretation indeed agrees with our equation 
(179) which can be written as 



0= / dyU{x,y) 



- Jdz U-\y,z)J{zyj+gc{y,y) 



(182) 
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That is, for a given J{x) ^ 0, one will solve as before the Hartree equation but with 
■'-'ion(x) — >■ fion(x) + J{x). Ouce the Hartree particle density uh is obtained, one obtains 
ikpc = uh + U~^oJ. Since Uo(iipc) always appears as a unit inside the Green's function 
Qc{x,x'), for convenience, we define 

Jc=U o{i(fic) — J + U oUh ■ 

Once (pc{x) is obtained, we shift by ipc and re-express the exponent in the integrand as 

—I[4>\ — i Jo0 =^ —I['^c + 0] — iJ°{<fc + (t>) 1 
and then expand around ipc- To do the expansion, we first rewrite (20) as 

G-^\x, x') = g-\x, x') + i h{x) 5{x - x') = g-\x, x') + V{x, x') , (183) 
where b = Uo(j). We may then write down 

G~^' = [I + ^coV] , 

and 

In {G-/) = In (^-^) + [OcoV]' . (184) 



k=l 

Note that 



Consequently, 

Trln(G'^i) = Trln(^;-^)+ J dxigc{xi,Xi){ib{xi)) 



~ j dxidx2Qc{xi.,X2)Qc{x2-,Xi){ih{xi)){ib{x2)) 

+ ^ j dxi. . .dxkQc{xk,Xi) . . .gc{xk-i,Xk){ih{xi)) . . .{ih{xk)) . (185) 

Note that in the final expression of the exponent of the integrand in (23) the terms linear 
in (or h) cancel out as one may verify and we arrive at 

-/[</?c + 0] - i Jo((^c + 0) = -]^^c°U°^c- Vc°b -]^hoU~'^oh+i:t\n{G~^^) 

-|--Trln(C/) — iJoifc — iJ°(t> 

= ^IV ln(C/) - \p^o Uoif^ + Tr hi{g-^) - iJoip, 



^ oo 

--boV;Kb + J2l^''Hvc]obi...obk , (186) 



k=3 
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where 



and 



I^''\cpc\obi...obk = 



i-1 



dxi... dxkQc{xk, xi)... Gc{xk-i, Xk) [{ib{^i)) ■ ■ ■ iib{^k))] ■ 

(187) 



We therefore have (based on the hnked-cluster theorem) 



PW^J] = -(pcoUoif^-Tr\n{g-^)+iJo(p, 



1 CX) 

+iTHn(p,W)-E^( 

n=l 



fe=3 



(188) 



where the subscript "conn." is used to indicate connected Feynman diagrams. 
Consequently, the effective action, defined by 

r[n] = ^W^J] - ^ Jo U'hj - Jon = ^W^[J] + -Jo U-^oJ - iJoip 

and with (/? = (/7c + (0) = </7c + <^ as well as with J — J^ — U oUh, can now be written as 

r[n] = -Trln {0'^) - J^oUh + ^uhoUoUh + ^Trln (v-^oU^ 



oo ^ 

-iJo(^- V— ( 

n=l 



2 2 

oo 

J2l^'\ip,]oh...ob, 



k=3 



(189) 



This expression should be contrasted with Eq. (45) where the true particle density is intro- 
duced as the natural variable. 

Prom the definition of (p (see Eq. (26) we have 

f D(j)(i(j)(x))e-^^^^-'^°^ /■ , 1. J Db Ub(y)) 6-^^^^+^^-^^°^^^+'!'^ 

M^) - ^ r ^Z\-L-.jo. - ^M^) + J dyU-' {x, y) ^ ^ 



which means that 



(190) 



(191) 
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On the basis of a simple replica argument^^, one sees that the above expression may be 
written as 



°° 1 /■ 

i(fi{x) = J^ — / dyU-\x,y){ 

oo „ 



k=3 

oo 



.k=3 



ih{y)). 



(192) 



where the n = term vanishes because / Dbe-2^°'^- °Hb{x) = 0. Eq. (192) implies that 
iip{x) at each point x is a functional of ipc- Returning to (189), this implies that one may 
express the effective action F in terms of iipc although its canonical argument is supposed 
to be n. The relation between n and i(fic is obtained through 



n — life + i(fi[i(pc] — U ^oJ[iipc] — uh + i^[i<fic 



(193) 



with iip[x) expressed as a functional of i(pc- This implies that given an iipc, we may obtain n 
through (i) the diagrammatic expansion of (192) to produce the corresponding iip, and (ii) 
equation (179) to find —U~^oJ. Since in general, we have 

dx (i(^(x)) ^ , 

one cannot expect the integral of the Hartree electron density to be Ng, the total number of 
electrons. Instead, we have 



Ng — J (ixn(x) = y (ix [nni'^) + i(p{x) 



(194) 



and J d:x.nH{^) is not necessarily an integer. In particular, when J d^KUni^) deviates sig- 
nificantly from Nc or J rfx |i(^(x)| ^ 1, forcing J d'^^nni^) = may lead to the occurrence 
of a self-consistent solution that differs significantly from the true solution. 

Note that nif(x) = —Qc{x,x) and in the absence of the source term nij(x) = i</7c(x). 
Furthermore, since {U^^oJ)^ oc VxJ(x), J njj(x) dx = J (i(pc{'x)) dx. The constraint (194) 
for the total number of electrons can also be written as 



Ne — J dxn(x) — J dx [njj(x) -|- i<^(x)] — J dx [i</?c(x) -|- i<^(x)] 



(195) 



The Hartree-like Green's function Qc{x,x') shown in (180) may be viewed as a func- 
tional of i(fic{x). When expressing Qc{x,x') using only single particle orbitals, we define 
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in Eq. (88) v(x) = fion(x) — A* + U-{i(fic):x.- For the evaluation of Qc{x,y), we solve first 
the eigensystem (93). The single-particle wave functions (91) associated with h are to 
be obtained self-consistently. Basically, one starts with a guess for the electronic den- 
sity iv'c(x) satisfying J (ix(i(/9c(x)) ^ N^., where is the number of electrons. One 
then obtains the single-particle wave functions, and then computes the corresponding 
Green's function Q^i obtains J cix(i(^(x)) and tunes the chemical potential to ensure that 
— J dx ^c(x, x) — Ng — J dx{i(p{x)). One then takes — ^c(x, x) in place of i</7c(x) in the next 
round of iteration until convergence is reached. 

The procedure of the saddle-point method is now obvious. One starts with an external 
potential, and then determines the Hartree electron density via n//(x) = —Qc{x,x), hh = 
i(Pc — U~^oJ and Eq. (195). This self- consistent procedure will also provide the physical 
electron density n(x). The ground state energy is obtained by using (189) to calculate the 
effective action, which is the ground state energy times /3 in the limit T — > 0. When one 
wishes to obtain the density functional r[n] at an electron density other than ut, one adds 
the source term into the potential fion(x) and then solves for the Hartree density, shown in 
(181), as outlined above. 

B. Remark on the single-particle limit 

The single electron limit for the Hartree-method is more complicated than for the method 
presented in section IV C. Because / nni'x.) dx is not necessarily 1 in the single electron hmit, 
due to Eq. (195), in general one needs to obtain rii self-consistently. The other issue is that 
the diagrammatic expansion contained in (189) does not cover all the diagrams for a given 
order of U. This means that one can't use the Hugenholtz diagram argument to eliminate 
vertex matrix element even when rii < 1. To see this point explicitly, we rewrite (189) as 

r[n] = -Tr In {Q-') + ( J - J,)onH + ^rino UoUh + ^Tr In (v;K [/) 



oo 

^ n! 

n=l 



^/(*^)[(^e]o6l...o6fe 



,fc=3 



Jon . (196) 
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Evidently, except for the last term above, the rest of the terms must constitute /3W^[J]. Since 
J — Jc = —UoUh, one can also write as 

fiW[J] = -Tr In (g;^) - Uho UoUh + ^uhoUoUh + ^Tr In (v'K 

I 

)conn. • 



-T- 

n=l 



5^/('=)[y.do6i...o6, 

.fc=3 

Diagrammatic expansion of the last two terms shows that the following diagrams 








of order f/^ are absent, when compared to the regular field theoretic perturbation calculation. 
This is not a disadvantage of the method. Instead, what our derivation shows is that the 
missing diagrams eventually will be compensated by the —nno UoriH term. However, it is 
obvious that the saddle-point formalism makes the single-electron limit hard to analyze. 

When ni > but ni ^ 1 at the low temperature limit, we know that Dc{x,y) = 
gc{x,y)Qc{y,x) will be of order rii. This is because 



Qc{x,y) = 5^0.(x)C(y)e-(^"-'^)(--^-) x 



-Ua) if < r, 



y 



[1 - ria) if > Ty 



and whenever < Ty, the propagator is of order rii. Since Dc{x,y) = Qc{x,y)Qc{y,x), 
one of the propagators in the product must be of order rii. In principle, one needs to 
solve for the occupation number rii of the lowest energy state using iipd'^) — iP^^oJ)^ = 
—Qc{x,x) and Eq. (195). Nevertheless, the correct one particle limit can be seen if one 
starts with a chemical potential /i such that ni ~ 0. In this case, we have at the J = limit 
i(Pc{x) = rini^) ~ as well as Dc{x,y) oc ni ^ 0. This way, the higher order exchange- 
correlation terms may be viewed as the having smaller contributions and one may control 
the accuracy by controlling the number of higher order terms included. Of course, one then 
has n(x) ^ «<^(x) and the condition J (ix(i(^(x)) = 1 — rii must be satisfied. 

C. Obtaining excitations using the Hartree method 

In general the excitations are determined by Eqs. (163) and (164). Under the Hartree 
formalism described in this section, the natural variable used is the Hartree density uh 
rather than the true particle density n^. One thus must transform the variable used in 
Eqs. (163) and (164) from ut to uh- We describe below how this can be achieved. 
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At the physical condition (J = 0), one has 



6n{x) 



dxidx2 



5V 



SnH(x2] 



SnH{x2) 


SJc{xi) 




n=nT ^Jc{xi) 







Note that Suh/SJc contains no zero mode because uh — SWh/SJc, where H/ij[Jc] = 
— Trln(^c~^), and S'^Wh/SJcSJc is known to be strictly negative from a theorem proved 
by Valiev and Fernando. ^° The strict negative-definiteness of 6'^W/SJSJ can also be inter- 
preted as the stability condition 



/ 



dx Sn{x) 5J{x) < , 



which means raising the local one-particle potential leads to an average decrease of the 
local particle concentration and vice versa. The strict negative-definiteness means that 
Suh/SJc contains no zero modes and is invertible. Also, because n — uh + i'^ci 5n/5Jc — 
Suh/^Jc + 5{i(p)/5Jc exists via diagrammatic expansion of i(p in terms of Jc- The existence 
of 5n/5Jc implies that 5Jc/5n is invertible (i.e., has no zero eigenvalues). Therefore, after 
multiplying the inverse of 5Jc/5n and Suh/SJc on both sides of the above equation, one has 

5V 



Snnix) 



0, 



n=nT 



which then leads to 



r^'\x,y) 



5n{x)5n{y) 



n=nT 



J 



dxidx2dyidy; 



SJc{xi) 6nH{x2) 



5^T 



5n{x) 5Jc{xi) SnH{x2)5nH{y2) ^Jcivi) 5n{y) 



n=nT 



Using Eq. (159), one may write r(^)(x, y, w) as the analytic continuation of the integral of 
r'^'^\x,y) over time. To achieve this goal, one first writes down 

r(2)(x,y,zi/„) = f d{Ty-T,)e'^-^'y-^^^V^^\x,y) 
Jo 

/•/3 

Jo 

= J dxidx2dyidy2f^\-iUn)g^l{-iJ^n)i'^'^^ (x2, y2, ii^n)gy%iJ^n)fy\iJ^n) 

= jdl^2dy2K\-lPn)T^^\^2,y2,ll^n)hf{lPn), (197) 
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where 



r(^)(x,y,ii/„) = / d{Ty-T, 
Jo 



5nH(x)5nH(y) 



hl\wn) = j dyi g^l{iun)fp{iun) . 
Therefore, Eq. (166) for finding the excitations becomes 



n=nT 



n=nT 



Jdzt^^\x,z,-uj)A{z)=0, (198) 



with A(z) given by 



A(z) = J dyhl{-uj)A{y). 

One therefore obtains the excitation energy via solving Eq. (198). Since hy{—uj) is invertible, 
one may also obtain A(z) via A(z) if desired. Within this framework, the first two terms 
of (189) correspond to the non-interacting part Fq of the effective action in section V. The 
protocol for obtaining the excitation energy is then the same as described in section V with 
the one exception that the ground state of the Hartree system is not the same as the ground 
state of the real system. 

VII. DISCUSSION AND FUTURE DIRECTIONS 

In this paper we focus on the auxiliary field method applied to the development of the 
density functional. It is natural to inquire into the physical meaning of the auxiliary, bosonic 
field introduced in eqs (15-20). In the path integral treatment of relativistic quantum 
electrodynamics, if one were to integrate out the photon field, one generates the current- 
current interaction which is quartic in fermionic field. When viewing this process backwards, 
one sees that the quartic fermionic interaction is disentangled by introducing the photon field. 
The field here thus plays a similar role to the photon field as it disentangles the quartic 
fermionic interaction term. As we will argue below, the variable is closely related to the 
"time" component of the photon field, i.e., the electric potential. 

To make the connection between the photon field and 0, let us first seek the nonrelativistic 
limit of the Lagrangian density, Cem — —F^''Fn^/{167r) , associated with the relativistic 
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photon field. Note that in the limit when the speed of light approaches infinity (removal of 
terms involving time derivative of the three- vector potential) , the Lagrangian density turns 
into (VAo)^/87r, which contains no quadratic derivative with respect to time. Therefore, at 
finite temperature (with it — )■ r) this implies that the exponent associated with the photon 
field path integral will behave as 



/ 



dt dx Cem ^ I dr I dx^^^^ ^ - f dr f dx^ AqV^Aq . (199) 



Setting C/(x — y) = e^/|x — y| [thus U ^{x,y) — —j^W'^5{-x — y)] and comparing (199) 
with Eqs. (15-16), we make the identification {iUo(j))/e = Aq, because now 

~4>oUo(t> ^Aq o{e\-^) oAo^- J dr j dx^ AqV^Aq . 

In the non-relativistic limit, the electric part and the magnetic part of electromagnetism 
are decoupled. Starting with a non-relativistic many-electron system, one may ask what is 
the quantum mechanical analog of the Poisson equation that forms the basis of electrostatics. 
It turns out that this is easily obtained via computing 6Z[J]/6J\ j^q in two different ways. 
First, upon taking the derivative of (15) and using Eqs. (16) and (17), one obtains 

6Z[J] 

j=o 



= -/3z,.o(^nx)^(x)) = _^^,^„MWM 



5J{x.) 

Second, while taking the derivative of (15), if using Eqs. (19) and (20) one arrives at 



SZ[J] 



= -pZj=om = -l3Zj=oe U-K{Ao) = /3Zj=ot^V^(Ao(x)) . 
,/=o ^"^^ 



5J(x) 

We therefore obtain the thermal quantum-mechanical analog of the Poisson equation 

V^(Ao(x)) = -47r(e^t(x)^(x)) ^ 

a result also obtained in reference 23. This connection to classical electrostatics is essen- 
tial since it provides the quantum-mechanical correspondence of an important ingredient 
in (bio) molecular interactions that have been extensively studied in the presence of di- 
electrics. ^^"^^ 

The UDF described in this paper is systematically constructed, uniquely determined, and 
in principle exact. However, in terms of real computations, one can only keep P^ terms up 
to some order in A. A natural question thus arises. How well will the truncated version 
work? In general, this question can only be answered with numerical results. However, by 



79 



providing theoretical arguments and comparisons to other approaches, we wish to convey 
that this method is hkely to produce good results and thus to attract computational efforts 
towards using the proposed approach. 

It is worth pointing out the relation between the expression (46) and the scheme^^ moti- 
vated by the renormalization- group. The vertex functions /(-^-'^ in (46) will contribute to the 
so-called Wocal approximation of reference 21. The absence of /^^^ manifests the absence of 
the correction term due to the bi- local contribution, as shown in reference 21. Furthermore, 
with the bifocal approximation^^ included, Polonyi and Sailer obtained an approximate en- 
ergy functional which corresponds exactly to our Fq -|- Fi . To reach an equivalent form of the 
proposed /-local approximation of reference 21, we simply keep terms up to F;/2- Therefore, 
our formulation provides an explicit means for achieving an /-local approximation without 
resorting to the Hellmann-Feynman theorem. 

As shown in (87), the propagator Vq can be expanded as 

Vo^U + UoDoo U + UoDqo UoDqo U+ ... , 

when U can be viewed as a small quantity and treated pcrturbatively. When one is not 
allowed to treat U as a small parameter (say in the strong coupling regime), or when one 
needs to treat as a small parameter instead, the conventional perturbative expansion 
in breaks down completely while our approach is still applicable. In the case when 
must be treated as small, we expand Vq as 

Vo = -D-^ - D-h U-^oD-^ - D-^o U-^oD-^o U'^oD'^ - ... . (200) 

And in this case, C/ ^ 1, our effective action expansion does have the Hartree term ^noUon 
as the leading order, followed by terms of order and then the expansion of T>o provides 
series in powers of U~^. Note that in this case, the exchange-correlation functional is not 
led by order U at all, but is led by order U^. This feature is not present in the conventional 
perturbative approach using U (or e^) as the expansion parameter. 

As mentioned earlier, there also exist different functional methods for many electron 
systems. For example, the exchange correlation functional outlined by Sham^^ is founded 
on the perturbative functional approach developed by Luttinger and Ward^^ or equiva- 
lently by Klein. A succinct review of the Luttinger- Ward/ Klein functional and its ap- 
plications can be found in reference 69. The Luttinger- Ward/ Klein functional yields the 
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grand-potential/ground-state-energy only when the functional argument is equal to the 
fully-interacting, physical, one-particle Green's function. Instead of allowing the physical, 
full, one-particle Green's function, Dahlen et al.^^ proposed to find that stationary point 
of the Luttinger- Ward/Klein functional while restricting the argument to the Hartree-Fock 
Green's functions or the Green's functions of non-interacting systems. However, even if the 
Luttinger- Ward/Klein functional is computed to all orders, the error in the value of the 
grand potential (or the ground state energy) due to restriction on the Green's functions 
remains unknown. 

Furthermore, it should be noted that when the functional argument is equal to the physi- 
cal onc-particlc Green's function, the Luttinger- Ward/Klein functional reaches a stationary 
point, not the minimum. This means that it is possible for the Klcin/Luttinger-Ward func- 
tional to assume an even lower value than the ground state energy (or the grand potential) 
when the functional argument deviates from the physical Green's function. In other words, 
the Klein/Luttinger-Ward functional only retains its meaning as the ground state energy (or 
grand potential) when the Green's function takes the value of the true (physical) Green's 
function. Our effective action expression of the energy functional, on the other hand, truly 
represents energy of the system. Our effective action energy functional, when no trunca- 
tion on the series is made, reaches its minimum when the electron density assumes the true 
(physical) density, and for any other v-representable density profile prescribed, it represents 
the lowest energy possible associated with that prescribed density profile. 

The method of Hedin^^ is largely identical to that of Luttinger and Ward.^^ This includes 
the fact that the energy functional reaches a stationary point, rather than the minimum, 
when the functional argument is the fully-interacting one-particle Green's function. How- 
ever, Hedin aims to replace the e^-based (bare Coulomb) perturbative expansion of the 
electron self energy by another expansion using a screened interaction W. Hedin expresses 
the electron self energy and the screened interaction as a functional of the electron Green's 
function of the interacting system. Interestingly, the first order result, also termed the GW 
approximation, of Hedin^^ has been shown^'^ to produce good results when compared to 
other density functional. This suggests that not treating as small might have some ad- 
vantage. It is worth pointing out that the first order term, Fi = — |Tr ln(P(^^o [/), of the 
UDF described here is equivalent to the celebrated GW approximation. 

Like the W propagator of Hedin, the T>o propagator introduced here also corresponds 
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to that of a screened interaction (see section IV D), thereby avoiding any possible infrared 
divergence associated with perturbative expansion based on bare Coulomb interactions. 
However, the screening associated with Vq is from the KS particles and thus keeps the same 
form no matter how many orders one wishes to include. This is different from that of Hedin's 
where the expression of W in terms of the electron Green's function changes with the order 
included. The other difference between the proposed approach and reference 31 is that the 
UDF proposed here depends on Jq, a function of three spatial variables (and possibly with 
one additional time variable), while the method of reference 31 expresses via W the electron 
self energy as a functional of the Green's function, a function of six spatial variables (and 
possibly with two additional time variables). 

It is well known that a loopwise expansion may also be viewed as an h expansion 
that is, an expansion of quantum-mechanical effects. By first integrating out the fermionic 
degrees of freedom completely, the proposed method is an expansion of bosonic loops formed 
by T>o propagators associated with the auxihary field b. The b field describes the potential 
produced by electron density fluctuations around Ug. Since the ground state charge density 
Ug captures the full quantum information of the ground state thanks to the HK theorem, 
one anticipates a weaker quantum effect associated with the auxiliary b field than with the 
fermionic field. This makes the auxiliary b field a suitable candidate for loop (or quantum 
effect) expansion, the approach pursued in this paper. 

Finally, let us remark on the issue of convexity. The full r[n] is supposed to be convex,^° 
thus guaranteeing a unique solution without any local minima when searching for the mini- 
mum of r[n]. However, in real computations only a finite number of terms of the effective 
action can be kept. This approximate/truncated expression may not warrant convexity and 
thus it is not guaranteed to be free of local minima while numerically searching for the 
ground state density Ug (or thermal averaged density ut at finite temperature) . In the near 
future, we plan to implement numerically the methods presented in this paper, and will 
describe in a separate publication the results obtained as well as the investigation on the 
issue of convexity. 
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